1 added + 32 removed + 44 modified, total 77 files
projects/slic/trunk
--- projects/slic/trunk/CMakeLists.txt 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/CMakeLists.txt 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,8 +1,9 @@
-#
-# CMake configuration for SLIC application.
-#
-# @author Jeremy McCormick <[log in to unmask]>
-#
+##########################################################
+# CMake configuration for SLIC application. #
+# #
+# @author Jeremy McCormick <[log in to unmask]> #
+# #
+##########################################################
# minimum cmake version
CMAKE_MINIMUM_REQUIRED( VERSION 2.8 )
@@ -11,9 +12,9 @@
PROJECT( SLIC )
# project version
-SET( SLIC_VERSION_MAJOR 3 )
-SET( SLIC_VERSION_MINOR 1 )
-SET( SLIC_VERSION_PATCH 7 )
+SET( SLIC_VERSION_MAJOR 4 )
+SET( SLIC_VERSION_MINOR 0 )
+SET( SLIC_VERSION_PATCH 0 )
# option for enabling and disabling the logging system
OPTION( ENABLE_SLIC_LOG "enable slic logging system" ON )
@@ -79,15 +80,16 @@
# find LCDD
FIND_PACKAGE( HEPPDT REQUIRED )
-IF ( NOT LCDD_FOUND )
+IF ( NOT HEPPDT_FOUND )
MESSAGE( FATAL_ERROR "HEPPDT was not found!" )
ELSE()
- MESSAGE( STATUS "Found HEPPDT: ${HEPPDT_ROOT}" )
+ MESSAGE( STATUS "Found HEPPDT: ${HEPPDT_DIR}" )
ENDIF()
# configure files
CONFIGURE_FILE( ${PROJECT_SOURCE_DIR}/include/PackageInfo.hh.in ${PROJECT_SOURCE_DIR}/include/PackageInfo.hh)
CONFIGURE_FILE( ${PROJECT_SOURCE_DIR}/include/Geant4VersionInfo.hh.in ${PROJECT_SOURCE_DIR}/include/Geant4VersionInfo.hh )
+CONFIGURE_FILE( ${PROJECT_SOURCE_DIR}/scripts/create_slic_tarball.sh.in ${PROJECT_SOURCE_DIR}/scripts/create_slic_tarball.sh )
# additional compilation flags for Geant4
ADD_DEFINITIONS( ${Geant4_DEFINITIONS} )
@@ -98,12 +100,9 @@
ADD_DEFINITIONS( -DSLIC_LOG )
ENDIF()
-# enable gprof when using debug build
-IF ( CMAKE_BUILD_TYPE MATCHES DEBUG )
- MESSAGE( STATUS "Enabling gprof in DEBUG build" )
- ADD_DEFINITIONS( "-pg" )
- SET( CMAKE_EXE_LINKER_FLAGS_DEBUG "${CMAKE_EXE_LINKER_FLAGS_DEBUG} -pg" )
-ENDIF()
+# debug flags
+SET ( CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS} -DDEBUG=1 -pg -ggdb" )
+SET ( CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS} -O3" )
# local includes
SET( SLIC_INCLUDE_DIR ${PROJECT_SOURCE_DIR}/include )
@@ -159,10 +158,10 @@
MESSAGE( STATUS "INSTALL_DOC is set to OFF. Documentation will not be created." )
ENDIF()
-# print CMake environment if enabled
-IF ( PRINT_CMAKE_VARIABLES STREQUAL "ON" )
- GET_CMAKE_PROPERTY( _variableNames VARIABLES )
- FOREACH ( _variableName ${_variableNames} )
- MESSAGE(STATUS "${_variableName}=${${_variableName}}" )
- ENDFOREACH()
-ENDIF()
+# install
+INSTALL( TARGETS slic slic DESTINATION bin )
+
+#get_cmake_property(_variableNames VARIABLES)
+#foreach (_variableName ${_variableNames})
+# message(STATUS "${_variableName}=${${_variableName}}")
+#endforeach()
projects/slic/trunk/cmake
--- projects/slic/trunk/cmake/FindHEPPDT.cmake 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/cmake/FindHEPPDT.cmake 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,7 +1,7 @@
FIND_PATH( HEPPDT_INCLUDE_DIR HepPDT/Version.hh ${HEPPDT_DIR}/include )
-FIND_LIBRARY( HEPPDT_LIBRARY HepPDT ${HEPPDT_DIR}/src/HepPDT/.libs/ )
-FIND_LIBRARY( HEPPID_LIBRARY HepPID ${HEPPDT_DIR}/src/HepPID/.libs/ )
+FIND_LIBRARY( HEPPDT_LIBRARY HepPDT ${HEPPDT_DIR}/lib/ )
+FIND_LIBRARY( HEPPID_LIBRARY HepPID ${HEPPDT_DIR}/lib/ )
SET( HEPPDT_LIBRARIES ${HEPPDT_LIBRARY} ${HEPPID_LIBRARY} )
projects/slic/trunk/include
--- projects/slic/trunk/include/CommandLineOption.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/CommandLineOption.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,7 +1,7 @@
// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/CommandLineOption.hh,v 1.3 2013-11-06 00:23:35 jeremy Exp $
-#ifndef SLIC_COMMANDLINEOPTION_HH
-#define SLIC_COMMANDLINEOPTION_HH
+#ifndef slic_CommandLineOption_hh_
+#define slic_CommandLineOption_hh_ 1
// STL
#include <map>
projects/slic/trunk/include
--- projects/slic/trunk/include/CommandLineProcessor.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/CommandLineProcessor.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,6 +1,6 @@
// $id: $
-#ifndef SLIC_COMMANDLINEPROCESSOR_HH
-#define SLIC_COMMANDLINEPROCESSOR_HH
+#ifndef slic_CommandLineProcessor_hh_
+#define slic_CommandLineProcessor_hh_ 1
// STL
#include <string>
@@ -88,8 +88,8 @@
* Print the version information.
*/
void printVersion() {
- log() << LOG::okay << PackageInfo::getFullApplicationString() << LOG::done;
- log() << LOG::okay << "Geant4 " << Geant4VersionInfo::getFullVersion() << LOG::done;
+ log() << LOG::okay << PackageInfo::getShortName() << " " << PackageInfo::getVersion() << LOG::done;
+ log() << LOG::okay << "Geant4 " << Geant4VersionInfo::getVersion() << LOG::done;
}
/**
projects/slic/trunk/include
--- projects/slic/trunk/include/CommandQueue.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/CommandQueue.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,6 +1,6 @@
// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/CommandQueue.hh,v 1.4 2013-11-06 00:23:35 jeremy Exp $
-#ifndef SLIC_COMMANDQUEUE_HH
-#define SLIC_COMMANDQUEUE_HH
+#ifndef slic_CommandQueue_hh_
+#define slic_CommandQueue_hh_ 1
// SLIC
#include "Module.hh"
projects/slic/trunk/include
--- projects/slic/trunk/include/EventDebugger.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/EventDebugger.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,129 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/EventDebugger.hh,v 1.7 2013-11-06 00:23:35 jeremy Exp $
-
-#ifndef SLIC_EVENTDEBUGGER_HH
-#define SLIC_EVENTDEBUGGER_HH 1
-
-// slic
-#include "Module.hh"
-#include "Singleton.hh"
-
-// geant4
-#include "globals.hh"
-
-// std
-#include <string>
-#include <vector>
-
-class G4Run;
-class G4Event;
-
-namespace slic {
-class EventDebuggerMessenger;
-
-/**
- * @class EventDebugger
- * @brief Manages debugging of individual G4Events.
- */
-class EventDebugger: public Module, public Singleton<EventDebugger> {
-public:
-
- typedef std::vector<G4int> DebugEventList;
-
-public:
-
- /**
- * Class constructor.
- */
- EventDebugger();
-
- /**
- * Class destructor.
- */
- virtual ~EventDebugger();
-
-public:
-
- /**
- * Turn event debugging on or off.
- * @param[in] e The event debugging setting.
- */
- void enableDebug(bool e = true);
-
- /**
- * Add an event number that should be debugged.
- * @param[in] eventNumber The event number to debug.
- */
- void addDebugEvent(G4int eventNumber);
-
- /**
- * Clear the event number list.
- */
- void clearDebugEvents();
-
- /**
- * Check if a given event should be debugged.
- * @param[in] eventNumber True if event should be debugged; false if not.
- */
- bool haveDebugEvent(G4int eventNumber) const;
-
- /**
- * Set a Geant4 macro to be executed on debug events.
- * @param[in] macro The Geant4 macro string to be executed on debug events.
- */
- void setDebugMacro(std::string macro);
-
- /**
- * Set the macro for cleaning up after debugging an event.
- * @param[in] macro The Geant4 macro string to be executed on debug events.
- */
- void setCleanupMacro(std::string macro);
-
- /**
- * Begin of run hook.
- * @param[in] run The Geant4 run.
- */
- void beginRun(const G4Run* run);
-
- /**
- * Begin of event hook.
- * @param[in] event The Geant4 event.
- */
- void beginEvent(const G4Event* event);
-
- /**
- * End of event hook.
- * @param[in] event The Geant4 event.
- */
- void endEvent(const G4Event* event);
-
-private:
-
- /**
- * Sort the event number list in ascending order.
- */
- void sortDebugEvents();
-
- /**
- * Execute the debug macro.
- */
- void execDebugMacro();
-
- /**
- * Execute the cleanup macro.
- */
- void execCleanupMacro();
-
-private:
-
- EventDebuggerMessenger* m_messenger;
- std::string m_debugMacro;
- std::string m_cleanupMacro;
- DebugEventList m_events;
- bool m_debugging;
- bool m_haveDebugMacro;
- bool m_haveCleanupMacro;
- bool m_forceDebugMode;
-};
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/EventDebuggerMessenger.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/EventDebuggerMessenger.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,60 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/EventDebuggerMessenger.hh,v 1.6 2013-11-06 00:23:35 jeremy Exp $
-
-#ifndef SLIC_EVENTDEBUGGERMESSENGER_HH
-#define SLIC_EVENTDEBUGGERMESSENGER_HH 1
-
-// geant4
-#include "G4UImessenger.hh"
-#include "G4UIdirectory.hh"
-#include "G4UIcommand.hh"
-
-class G4UIcmdWithAString;
-class G4UIcmdWithABool;
-
-namespace slic {
-
-/**
- * @class EventDebuggerMessenger
- * @brief G4UImessenger for the EventDebugger.
- */
-class EventDebuggerMessenger: public G4UImessenger {
-
-public:
-
- /**
- * Class constructor.
- */
- EventDebuggerMessenger();
-
- /**
- * Class destructor.
- */
- virtual ~EventDebuggerMessenger();
-
-public:
-
- /**
- * Process the macro command.
- * @param[in] cmd The macro command.
- * @param[in] newVals The command arguments.
- */
- void SetNewValue(G4UIcommand* cmd, G4String newVals);
-
-private:
-
- /**
- * Define the Geant4 macro commands.
- */
- void defineCommands();
-
-private:
- G4UIdirectory* m_debugDir;
- G4UIcmdWithAString* m_debugMacroCmd;
- G4UIcmdWithAString* m_cleanupMacroCmd;
- G4UIcommand* m_clearDebugEventsCmd;
- G4UIcmdWithABool* m_enableDebugCmd;
- G4UIcommand* m_addDebugEventsCmd;
-};
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/EventSource.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/EventSource.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -45,7 +45,7 @@
/**
* Dump information about the current event.
*/
- virtual void dumpCurrentEvent() = 0;
+ virtual void printCurrentEvent() = 0;
/**
* Reset the event source.
projects/slic/trunk/include
--- projects/slic/trunk/include/EventSourceManager.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/EventSourceManager.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -19,6 +19,7 @@
#include "G4UnitsTable.hh"
namespace slic {
+
/**
* @class EventSourceManager
* @brief This singleton is responsible for managing physics event sources.
@@ -36,8 +37,14 @@
public:
+ /**
+ * Class destructor.
+ */
virtual ~EventSourceManager();
+ /**
+ * Class constructor.
+ */
EventSourceManager();
public:
@@ -274,61 +281,69 @@
*/
void setZSmearing(const G4double zSmear) {
m_zSmearingParam = zSmear;
- G4cout << "Set m_zSmearingParam = " << m_zSmearingParam << G4endl;}
+ G4cout << "Set m_zSmearingParam = " << m_zSmearingParam << G4endl;
+ }
- /**
- * Get the Z smearing parameter.
- * @return The Z smearing parameter.
- */
- G4double getZSmearing() {
- return m_zSmearingParam;
- }
+ /**
+ * Get the Z smearing parameter.
+ * @return The Z smearing parameter.
+ */
+ G4double getZSmearing() {
+ return m_zSmearingParam;
+ }
- public:
+ void enablePrintEvent(G4bool enablePrintEvent) {
+ _enablePrintEvent = enablePrintEvent;
+ }
- virtual void GeneratePrimaryVertex(G4Event* evt);
+public:
- private:
+ void generateNextEvent(G4Event* evt);
- // generator messenger
- GeneratorMessenger* m_messenger;
+private:
- // current event source
- EventSource* m_currentEventSource;
+ // generator messenger
+ GeneratorMessenger* m_messenger;
- // ParticleGun source. Always enabled.
- EventSource* m_particleGunSource;
+ // current event source
+ EventSource* m_currentEventSource;
- // filename
- std::string m_filename;
- bool m_fileIsSet;
- bool m_newFilename;
+ // ParticleGun source. Always enabled.
+ EventSource* m_particleGunSource;
- // num events generated on current generator
- int m_ngen;
+ // filename
+ std::string m_filename;
+ bool m_fileIsSet;
+ bool m_newFilename;
- // num events to skip
- unsigned int m_nskip;
+ // num events generated on current generator
+ int m_ngen;
- // need to setup a new source new run?
- bool m_newSource;
+ // num events to skip
+ unsigned int m_nskip;
- // 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;
+ // need to setup a new source new run?
+ bool m_newSource;
- // source type
- ESourceType m_sourceType;
+ // gen strings
+ static std::string STDHEP;
+ static std::string LCIO;
+ static std::string GPS;
+ static std::string PARTICLE_GUN;
+ static std::string UNKNOWN;
- // Lorentz transformation angle which will be applied to event.
- G4double m_lorentzTransformationAngle;
+ // source type
+ ESourceType m_sourceType;
- // Parameter which will be used to smear vertex Z position.
- G4double m_zSmearingParam;
- };
- }
+ // Lorentz transformation angle which will be applied to event.
+ G4double m_lorentzTransformationAngle;
+ // Parameter which will be used to smear vertex Z position.
+ G4double m_zSmearingParam;
+
+ G4bool _enablePrintEvent;
+
+};
+}
+
#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/EventSourceWithInputFile.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/EventSourceWithInputFile.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -95,7 +95,7 @@
/**
* Dump the current event to std::cout.
*/
- virtual void dumpCurrentEvent() = 0;
+ virtual void printCurrentEvent() = 0;
protected:
projects/slic/trunk/include
--- projects/slic/trunk/include/GPSEventSource.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/GPSEventSource.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -40,7 +40,7 @@
/**
* Dump the current event to std::cout.
*/
- void dumpCurrentEvent();
+ void printCurrentEvent();
/**
* Get the Geant4 general particle source class.
projects/slic/trunk/include
--- projects/slic/trunk/include/Geant4VersionInfo.hh.in 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/Geant4VersionInfo.hh.in 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,6 +1,10 @@
#ifndef Geant4VersionInfo_hh
#define Geant4VersionInfo_hh 1
+// LCDD
+#include "lcdd/util/StringUtil.hh"
+
+// STL
#include <string>
/**
@@ -9,23 +13,35 @@
*/
class Geant4VersionInfo
{
- private:
+private:
- /**
- * Class constructor.
- */
- Geant4VersionInfo() {;}
+ /**
+ * Class constructor.
+ */
+ Geant4VersionInfo() {;}
- public:
+public:
- /**
- * Get the Geant4 version string in format "##.##.##"
- * or "[release].[version].[patch]".
- * @return The Geant4 version string.
- */
- static std::string getFullVersion() {
- return "@Geant4_VERSION@";
- }
+ /**
+ * Get the Geant4 version string in format "##.##.##"
+ * or "[release].[version].[patch]".
+ * @return The Geant4 version string.
+ */
+ static std::string getVersion() {
+ return "@Geant4_VERSION@";
+ }
+
+ /**
+ * Get the Geant4 version string in format "##_##_##"
+ * or "[release]_[version]_[patch]".
+ * @return The Geant4 version string.
+ */
+ static std::string getVersion(std::string& sep, int width = 2) {
+ std::string version = StringUtil::pad(std::string("@Geant4_VERSION_MAJOR@"), width);
+ version += sep + StringUtil::pad(std::string("@Geant4_VERSION_MINOR@"), width);
+ version += sep + StringUtil::pad(std::string("@Geant4_VERSION_PATCH@"), width);
+ return version;
+ }
};
#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/GeneratorMessenger.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/GeneratorMessenger.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -56,7 +56,7 @@
G4UIcommand* m_resetCmd;
G4UIcmdWithAnInteger* m_skipEventsCmd;
G4UIcommand* m_filenameCmd;
- G4UIcommand* m_dumpCurrentEventCmd;
+ G4UIcommand* m_enablePrintEventCmd;
G4UIcommand* m_printNumEventsGeneratedCmd;
G4UIcmdWithADoubleAndUnit* m_setLorentzTransformationAngleCmd;
G4UIcmdWithADouble* m_setZSmearingParameterCmd;
projects/slic/trunk/include
--- projects/slic/trunk/include/LcioEventSource.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/LcioEventSource.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,92 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/LcioEventSource.hh,v 1.10 2013-11-06 00:23:35 jeremy Exp $
-
-#ifndef SLIC_LCIOEVENTSOURCE_HH
-#define SLIC_LCIOEVENTSOURCE_HH 1
-
-// slic
-#include "EventSourceWithInputFile.hh"
-#include "Module.hh"
-#include "LcioMcpFilter.hh"
-#include "LcioManager.hh"
-#include "LcioMcpManager.hh"
-
-// lcio
-#include "EVENT/LCEvent.h"
-#include "EVENT/LCCollection.h"
-#include "IOIMPL/LCFactory.h"
-#include "IO/LCReader.h"
-
-namespace slic {
-/**
- * @class LcioEventSource
- * @brief Implements event generation from the MCParticle collection of an LCIO file.
- */
-class LcioEventSource: public EventSourceWithInputFile {
-
-public:
-
- /**
- * Class constructor.
- * @param[in] fname The file name.
- */
- LcioEventSource(const std::string& filename = "");
-
- /**
- * Class destructor.
- */
- virtual ~LcioEventSource();
-
-public:
-
- /**
- * Open the current file.
- */
- void open();
-
- /**
- * Close the current file.
- */
- void close();
-
- /**
- * Read the next event.
- */
- void readNextEvent();
-
- /**
- * Dump the current event.
- */
- void dumpCurrentEvent();
-
- /**
- * Generate a Geant4 event.
- * @param[in] anEvent The target event.
- */
- void generate(G4Event* anEvent);
-
- /**
- * Begin run hook.
- * @param[in] aRun The run.
- */
- void beginRun(const G4Run* aRun);
-
- /**
- * Begin event hook.
- * @param[in] anEvent The event.
- */
- void beginEvent(const G4Event* anEvent);
-
- /**
- * Get the current MCParticle collection.
- * @return The current MCParticle collection.
- */
- EVENT::LCCollection* getCurrentMcpLCCollection();
-
-private:
- LcioMcpFilter* m_filter;
- IO::LCReader* m_reader;
- EVENT::LCEvent* m_event;
-};
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/LcioHitsCollectionBuilder.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/LcioHitsCollectionBuilder.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -193,9 +193,6 @@
// store momentum
bool m_storeMomentum;
-
- // ptr to Lcio Mcp manager for lkp of Mcp data
- LcioMcpManager* m_mcpManager;
};
}
projects/slic/trunk/include
--- projects/slic/trunk/include/LcioManager.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/LcioManager.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -9,11 +9,8 @@
#include "lcdd/hits/TrackerHit.hh"
// slic
-#include "LcioPrimaryGenerator.hh"
-#include "LcioMcpManager.hh"
#include "LcioHitsCollectionBuilder.hh"
#include "LcioFileNamer.hh"
-#include "Trajectory.hh"
#include "Singleton.hh"
// lcio
@@ -36,7 +33,6 @@
class G4Run;
class G4VHitsCollection;
-class G4TrajectoryContainer;
namespace slic {
@@ -99,18 +95,6 @@
void setAutonameFields(const std::vector<std::string>& fields);
/**
- * Get the output path for the LCIO files.
- * @return The output path for the LCIO files.
- */
- const std::string& getPath() const;
-
- /**
- * Get the output file name.
- * @return The output file path.
- */
- const std::string& getFilename() const;
-
- /**
* Get the path plus the filename and optionally include the extension.
* @param[in] withExtension Whether or not to include the file extension.
* @return The full output file path.
@@ -156,11 +140,6 @@
IMPL::LCEventImpl* createLCEvent();
/**
- * Create the final output particle collction.
- */
- void createFinalMcpCollection();
-
- /**
* Add a collection to the LCEvent.
* @param[in] lcsevent The LCEvent.
* @param[in] collection The collection.
@@ -175,11 +154,6 @@
void addCollection(EVENT::LCCollection* collection, const std::string& name);
/**
- * Add the initial particle collection.
- */
- void addInitialMCParticleCollection();
-
- /**
* End of event hook to write out an LCIO file.
* @param[in] The physics event.
*/
@@ -233,26 +207,10 @@
}
/**
- * Get the LCIO generator.
- * @return The LCIO generator.
- */
- inline LcioPrimaryGenerator* getGenerator() {
- return m_eventGenerator;
- }
-
- /**
- * Get the MCParticle manager.
- * @return The MCParticle manager.
- */
- inline LcioMcpManager* getMcpManager() {
- return m_mcpManager;
- }
-
- /**
* Get the LcioHitsCollectionBuilder for creating output hit collections.
* @return The LcioHitsCollectionBuilder.
*/
- inline LcioHitsCollectionBuilder* getHCBuilder() {
+ inline LcioHitsCollectionBuilder* getHitsCollectionBuilder() {
return m_HCBuilder;
}
@@ -272,6 +230,18 @@
m_enableDumpEvent = p;
}
+ /**
+ * Enable writing complete event to MCParticle collection.
+ * @param[in] writeCompleteEvent Set to true to write complete event.
+ */
+ void setWriteCompleteEvent(G4bool writeCompleteEvent) {
+ _writeCompleteEvent = writeCompleteEvent;
+ }
+
+ G4bool getWriteCompleteEvent() {
+ return _writeCompleteEvent;
+ }
+
private:
/**
@@ -320,16 +290,8 @@
*/
void addActiveSubdetectors();
- /**
- * Reset the instance variables of this class.
- */
- void reset();
-
private:
- // final Mcp collection with all info from McpInitial, Trajectories, Primaries
- IMPL::LCCollectionVec* m_McpFinalColl;
-
// writer
IO::LCWriter* m_writer;
@@ -340,11 +302,8 @@
std::string m_filename;
std::string m_path;
- // Mcp manager
- LcioMcpManager* m_mcpManager;
-
// event generator from MCP Coll
- LcioPrimaryGenerator* m_eventGenerator;
+ //LcioPrimaryGenerator* m_eventGenerator;
// creation of HCs
LcioHitsCollectionBuilder* m_HCBuilder;
@@ -369,7 +328,7 @@
bool m_enableDumpEvent;
bool m_writerIsOpen;
bool m_usingAutoname;
- //bool m_abortCurrentRun;
+ bool _writeCompleteEvent;
};
}
projects/slic/trunk/include
--- projects/slic/trunk/include/LcioMcpFactory.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/LcioMcpFactory.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,175 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/LcioMcpFactory.hh,v 1.11 2013-11-06 00:23:35 jeremy Exp $
-
-#ifndef SLIC_LCIOMCPFACTORY_HH
-#define SLIC_LCIOMCPFACTORY_HH 1
-
-// lcio
-#include "EVENT/LCCollection.h"
-#include "EVENT/LCEvent.h"
-#include "EVENT/MCParticle.h"
-#include "IMPL/MCParticleImpl.h"
-
-// slic
-#include "Module.hh"
-
-class G4PrimaryParticle;
-class G4TrajectoryContainer;
-class G4Event;
-
-namespace slic {
-class LcioMcpManager;
-class Trajectory;
-
-/**
- * @class LcioMcpFactory
- * @brief Creates LCIO MCParticles from saved Geant4/SLIC runtime information.
- */
-class LcioMcpFactory: public Module {
-
-public:
-
- /**
- * Class constructor.
- * @param[in] manager The LcioMcpManager used for particle utilities.
- */
- LcioMcpFactory(LcioMcpManager* manager);
-
- /**
- * Class destructor.
- */
- virtual ~LcioMcpFactory();
-
-public:
-
- /**
- * Create a new, empty MCParticle collection.
- * @param[in] collectionName The name of the collection.
- * @param[in] errorOnExists Whether to throw an error if the collection already exists (otherwise will overwrite).
- * @return The new, empty LCCollection.
- */
- EVENT::LCCollection* createEmptyMcpCollection(const std::string& collectionName, bool errorOnExist = true);
-
- /**
- * Create the final MCParticle output collection from the event.
- * @param[in] event The current Geant4 event.
- */
- void createFinalMcpCollection(const G4Event* event);
-
-private:
-
- /**
- * Create the final MCParticle vector for an event using the MCParticle vector
- * from the StdHep input.
- * @param[in] mcpVecInitial The initial LCIO MCParticle collection.
- */
- void createFinalCollectionFromInitial(EVENT::LCCollection* mcpVecInitial);
-
- /**
- * Create the output MCParticle collection using trajectory data only.
- * @param[in] trajectories The trajectory container.
- */
- void createFinalCollectionFromTrajectories(G4TrajectoryContainer* trajectories);
-
- /**
- * Create the MCParticle daughters of a trajectory that has the given track ID.
- * @param[in] parentParticle The parent MCParticle.
- * @param[in] parentID The parent's track ID.
- */
- void createDaughters(IMPL::MCParticleImpl* parentParticle, int parentID);
-
- /**
- * Create an MCParticle from a trajectory only.
- * @param[in] The trajectory.
- * @return The new MCParticle.
- */
- IMPL::MCParticleImpl* fromTrajectory(Trajectory* trajectory);
-
- /**
- * Create MCParticles recursively given an initial starting particle.
- * @param[in] The start point for creating daughter particles recursively.
- * @return The initial top particle with daughters added.
- */
- IMPL::MCParticleImpl* fromInitialRecurse(EVENT::MCParticle* particle);
-
- /**
- * Create an MCParticle from an initial MCParticle, only.
- * @param[in] mcpInit The initial MCParticle.
- * @return The new MCParticle.
- */
- IMPL::MCParticleImpl* fromInitialOnly(EVENT::MCParticle* mcpInit);
-
- /**
- * Create an MCParticle from a G4PrimaryParticle.
- * @param[in] The input G4PrimaryParticle.
- * @param[in] The corresponding initial MCParticle.
- * @return The new MCParticle.
- */
- IMPL::MCParticleImpl* fromPrimary(G4PrimaryParticle* primary, EVENT::MCParticle* particle);
-
- /**
- * Create an MCParticle from a G4PrimaryParticle and associated, initial MCParticle.
- * @param[in] primary The input G4PrimaryParticle.
- * @param[in] initial The associated initial MCParticle.
- * @return The new MCParticle.
- */
- IMPL::MCParticleImpl* fromInitialAndPrimary(G4PrimaryParticle* primary, EVENT::MCParticle* initial);
-
- /**
- * Create an MCParticle from a trajectory and an input MCParticle.
- * @param[in] trajectory The G4Trajectory.
- * @param[in] initial The initial, associated MCParticle.
- * @return The new MCParticle.
- */
- IMPL::MCParticleImpl* fromInitialAndTrajectory(Trajectory* trajectory, EVENT::MCParticle* initial);
-
- /**
- * Create an MCParticle from a G4PrimaryParticle via shallow copy, e.g. daughters are not created automatically.
- * @param[in] The G4PrimaryParticle.
- * @return The new MCParticle.
- */
- IMPL::MCParticleImpl* fromPrimaryShallowCopy(G4PrimaryParticle* primary);
-
- /**
- * Shallow copy an MCParticle.
- * @param[in] particle The MCParticle to copy.
- * @return The new MCParticle.
- */
- IMPL::MCParticleImpl* shallowCopy(EVENT::MCParticle* particle);
-
- /**
- * Create the daughter MCParticles from the G4PrimaryParticle and initial MCParticle
- * and add them to the parent.
- * @param[in] primary The starting G4PrimaryParticle.
- * @param[in] initial The associated, initial MCParticle.
- */
- void daughtersFromPrimary(G4PrimaryParticle* primary, EVENT::MCParticle* initial, IMPL::MCParticleImpl* mcpPar);
-
- /**
- * Add daughters to an MCParticle based on the associated initial MCParticle.
- * @param[in] newParticle The new MCParticle to be updated.
- * @param[in] initialParticle The initial MCParticle containing daughter information.
- */
- void addDaughtersFromInitial(IMPL::MCParticleImpl* newParticle, EVENT::MCParticle* initialParticle);
-
- /**
- * Create an LCGenericObject collection with MCParticle endpoint energies.
- * The order will match that of the MCParticle collection.
- * @param[in] particleCollection The MCParticle collection.
- */
- void fillEndPointEnergy(IMPL::LCCollectionVec* particleCollection);
-
- /**
- * Copy spin and color flow information from one particle to another.
- * @param[in] sourceParticle The input MCParticle.
- * @param[out] targetParticle The output MCParticle with updated spin and color flow information.
- */
- void copySpinAndColorFlow(EVENT::MCParticle* sourceParticle, EVENT::MCParticle* targetParticle);
-
-private:
- LcioMcpManager* m_manager;
- IMPL::LCCollectionVec* m_finalColl;
- G4TrajectoryContainer* m_currentTrajectoryContainer;
-};
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/LcioMcpFilter.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/LcioMcpFilter.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,69 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/LcioMcpFilter.hh,v 1.9 2013-11-06 00:23:35 jeremy Exp $
-
-#ifndef SLIC_LCIOMCPFILTER_HH
-#define SLIC_LCIOMCPFILTER_HH 1
-
-// slic
-#include "Module.hh"
-
-// lcio
-#include "EVENT/LCEvent.h"
-#include "EVENT/LCCollection.h"
-#include "EVENT/MCParticle.h"
-#include "IMPL/MCParticleImpl.h"
-
-namespace slic {
-/**
- * @class LcioMcpFilter
- * @brief Filters simulation particles from an MCParticle collection.
- */
-class LcioMcpFilter: public Module {
-
-public:
-
- /**
- * Class constructor.
- */
- LcioMcpFilter();
-
- /**
- * Class destructor.
- */
- virtual ~LcioMcpFilter();
-
-public:
-
- /**
- * Create a new MCParticle collection with only generator particles included.
- * @param[in] inputCollection The input MCParticle collection.
- * @return The filtered, output collection.
- */
- EVENT::LCCollection* filterSimParticles(EVENT::LCCollection* inputCollection);
-
- /**
- * Create a new MCParticle and add it to the target collection, unless it was already added.
- * @param[in] particle The input MCParticle to be copied.
- * @param[out] targetCollection The collection to be updated.
- * @return The new MCParticle.
- */
- IMPL::MCParticleImpl* makeParticle(EVENT::MCParticle* particle, EVENT::LCCollection* targetCollection);
-
- /**
- * Add simulation daughters from a parent particle to the target.
- * @param[out] targetParticle The target particle to update.
- * @param[in] parentParticle The source particle with the daughters.
- */
- void addSimDaughters(EVENT::MCParticle* targetParticle, EVENT::MCParticle* parentParticle, EVENT::LCCollection* mcpColl);
-
- /**
- * Check if a particle was already added to the list.
- * @return True if particle is already in list; false if not.
- */
- bool wasAdded(EVENT::MCParticle* particle);
-
-private:
- std::vector<EVENT::MCParticle*> m_mcpList;
-};
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/LcioMcpManager.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/LcioMcpManager.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,230 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/LcioMcpManager.hh,v 1.57 2013-11-06 00:23:35 jeremy Exp $
-
-#ifndef SLIC_LCIOMCPMANAGER_HH
-#define SLIC_LCIOMCPMANAGER_HH 1
-
-// slic
-#include "Module.hh"
-#include "LcioMcpMaps.hh"
-#include "LcioPrimaryGenerator.hh"
-#include "Singleton.hh"
-
-//#include "MCParticleMap.hh"
-
-// lcio
-#include "EVENT/LCCollection.h"
-#include "EVENT/LCEvent.h"
-#include "EVENT/MCParticle.h"
-#include "IMPL/MCParticleImpl.h"
-
-// geant4
-#include "globals.hh"
-#include "G4PrimaryParticle.hh"
-
-// std
-#include <string>
-
-class G4Event;
-class G4TrajectoryContainer;
-
-namespace slic {
-
-class LcioMcpMessenger;
-class LcioMcpFactory;
-class LcioMcpPrinter;
-class Trajectory;
-
-/**
- * @class LcioMcpManager
- * @brief This singleton manages creation of the LCIO MCParticle output collection from transient
- * G4Event and slic data.
- */
-class LcioMcpManager: public Module, public Singleton<LcioMcpManager> {
-
-public:
- // Mcp collections
- typedef std::map<std::string, EVENT::LCCollection*> McpCollectionMap;
-
-public:
- virtual ~LcioMcpManager();
- LcioMcpManager();
-
-public:
-
- /**
- * End of event hook.
- * @param[in] anEvent The event.
- */
- void endEvent(const G4Event* anEvent);
-
- /**
- * Get the MCParticle collection generated that was generated from an input StdHep event.
- * This will create an empty collection if it does not exist.
- * @return The MCParticle collection from the input StdHep event.
- */
- EVENT::LCCollection* getInitialMcpCollection();
-
- /**
- * Get the MCP collection generated at EndOfEvent.
- * Creates empty collection if doesn't exist.
- */
- EVENT::LCCollection* getFinalMcpCollection();
-
- /**
- * Set the minimum tracking distance for creating G4PrimaryParticles and G4PrimaryVertex objects
- * during event generation.
- * @param[in] minimumTrackingDistance The minimum tracking distance.
- */
- void setMinimumTrackingDistance(double minimumTrackingDistance);
-
- /**
- * Get the minimum tracking distance.
- * @return The minimum tracking distance.
- */
- double getMinimumTrackingDistance();
-
- /**
- * Reset the manager's state. To be called after done processing a single event.
- */
- void reset();
-
- /**
- * Set the initial MCParticle collection created from event generation.
- * @param[in] initialCollection The initial MCParticle collection.
- */
- void setInitialMcpCollection(EVENT::LCCollection* initialCollection);
-
- /**
- * Register by name an LCCollection containing MCParticles.
- * @param[in] name The name of the collection.
- * @param[in] collection The collection to register.
- */
- void registerMcpCollection(const std::string& name, EVENT::LCCollection* collection);
-
- /**
- * Find the first MCParticle collection in the LCevent.
- * @param[in] event The LCEvent to search.
- * @return The found collection or null if not found.
- */
- EVENT::LCCollection* findParticleCollection(EVENT::LCEvent* event);
-
- /*
- * Find a registered MCParticle collection by name.
- * @param[in] name The name of the collection to find.
- * @return The found collection or null if not found.
- */
- EVENT::LCCollection* findParticleCollection(const std::string& name);
-
- /**
- * Using the LcioMcpPrinter, print out the named collection from the current event.
- * @param[in] name The name of the collection to print out.
- */
- void printMcpCollection(const std::string& name);
-
- /**
- * Using the LcioMcpPrinter, print out the named collection from the current event.
- * @param[in] name The name of the collection to print out.
- */
- void printMcpCollection(const std::string& collName, EVENT::LCCollection* coll);
-
- /**
- * Get the MCParticle data maps.
- * @return The MCParticle maps.
- */
- LcioMcpMaps* getMaps();
-
- /**
- * Enable printing of the final MCParticle collection.
- * @param[in] p The enable setting.
- */
- void enablePrintFinal(bool p = true);
-
- /**
- * Enable printing of the initial MCParticle collection.
- * @param[in] p The enable setting.
- */
- void enablePrintInitial(bool p = true);
-
- /**
- * Begin of event hook.
- * @param[in] anEvent The event.
- */
- void beginEvent(const G4Event* anEvent);
-
- /**
- * Set whether to write out the initial MCParticle collection to the event
- * as a supplementary collection.
- * @param[in] b Set to true in order to write this collection; false to not write.
- */
- inline void setWriteInitialMCParticleCollection(bool b) {
- m_writeInitialMCParticleCollection = b;
- }
-
- /**
- * Get whether to write out the intial MCParticle collection to the event.
- * @return True if initial MCParticle collection will be written; false if not.
- */
- inline bool writeInitialMCParticleCollection() {
- return m_writeInitialMCParticleCollection;
- }
-
- /**
- * Get the index of a particle in a collection or -1 if not present.
- * @param[in] collection The LCCollection to search.
- * @param[in] particle The MCParticle to find.
- * @return The index of the MCParticle or -1 if not found.
- */
- int getParticleIndex(EVENT::LCCollection* collection, EVENT::MCParticle* particle);
-
-private:
-
- /**
- * Delete the current, initial MCParticle collection.
- */
- void deleteInitialMcpCollection();
-
- /*
- * Create an empty MCParticle collection indexed by name and register it in the map.
- * @param[in] name The name of the collection.
- * @return The new MCParticle collection.
- */
- EVENT::LCCollection* createParticleCollection(const std::string& name);
-
- // clear the maps
- void clearMaps();
-
-private:
-
- // current LCEvent
- IMPL::LCEventImpl* m_currentLCEvent;
-
- // printer for MCP debugging
- LcioMcpPrinter* m_printer;
-
- // user settings
- double m_vertexIsNotEndpointOfParentTolerance;
- double m_minimumTrackingDistance;
-
- // G4 messenger for Mcp handling
- LcioMcpMessenger* m_messenger;
-
- // map of Mcp colls being tracked by the manager
- McpCollectionMap m_mcpColls;
-
- /* Factory for creating MCParticles from Geant4 runtime data. */
- LcioMcpFactory* m_factory;
-
- /* Collection of maps with MCParticle information. */
- LcioMcpMaps* m_maps;
-
- bool m_enablePrintFinal;
- bool m_enablePrintInitial;
- bool m_writeInitialMCParticleCollection;
-
-public:
- static const double m_NAN;
- static const double DEFAULT_MIN_TRACKING_DISTANCE;
-};
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/LcioMcpMaps.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/LcioMcpMaps.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,170 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/LcioMcpMaps.hh,v 1.8 2013-11-06 00:23:35 jeremy Exp $
-
-#ifndef SLIC_LCIOMCPMAPS_HH
-#define SLIC_LCIOMCPMAPS_HH 1
-
-// slic
-#include "Module.hh"
-
-// lcio
-#include "IMPL/MCParticleImpl.h"
-#include "EVENT/MCParticle.h"
-
-// geant4
-#include "globals.hh"
-
-// stl
-#include <map>
-
-// geant4
-class G4PrimaryParticle;
-
-namespace slic {
-
-// slic
-class LcioMcpManager;
-class Trajectory;
-
-/**
- * @class LcioMcpMaps
- * @brief Data maps that relate MCParticles to various other runtime objects.
- */
-class LcioMcpMaps: public Module {
-
-public:
-
- /**
- * Class constructor.
- * @param[in] particleManager The MCParticle Manager.
- */
- LcioMcpMaps(LcioMcpManager* particleManager);
-
- /**
- * Class destructor.
- */
- virtual ~LcioMcpMaps();
-
-public:
-
- /**
- * Map of track IDs to MCParticle.
- */
- typedef std::map<G4int, IMPL::MCParticleImpl*> TrackToMcpMap;
-
- /**
- * Map of MCParticle to G4PrimaryParticle.
- */
- typedef std::map<EVENT::MCParticle*, G4PrimaryParticle*> McpToPrimaryMap;
-
- /**
- * Map of initial MCParticle from StdHep to final MCParticle.
- */
- typedef std::map<EVENT::MCParticle*, IMPL::MCParticleImpl*> InitMcpToMcpMap;
-
-public:
-
- /**
- * Map an initial MCParticle to a G4PrimaryParticle.
- * @param[in] particle The MCParticle.
- * @param[in] primary The G4PrimaryParticle.
- */
- void addParticleToPrimary(EVENT::MCParticle* particle, G4PrimaryParticle* primary);
-
- /**
- * Map a track ID to an MCParticle.
- * @param[in] trackID The track ID.
- * @param[in] particle The MCParticle.
- */
- void addTrackToMCParticle(G4int trackID, IMPL::MCParticleImpl* particle);
-
- /**
- * Map an initial MCParticle to a final output MCParticle.
- * @param[in] initial The initial MCParticle.
- * @param[in] final The final MCParticle.
- */
- void addInitialParticleToFinal(EVENT::MCParticle* initial, IMPL::MCParticleImpl* final);
-
- /**
- * Find an MCParticle from the associated trackID.
- * @param[in] trackID The track ID.
- * @return The MCParticle or null if not found.
- */
- IMPL::MCParticleImpl* findParticle(G4int trackID) const;
-
- /**
- * Find a trackID from a final MCParticle.
- * @param[in] particle The final MCParticle.
- * @return The track ID or -1 if not found.
- */
- G4int findTrack(EVENT::MCParticle* particle) const;
-
- /**
- * Find a G4PrimaryParticle from an MCParticle.
- * @param[in] particle The MCParticle.
- * @return The G4PrimaryParticle or null if not found.
- */
- G4PrimaryParticle* findPrimary(EVENT::MCParticle* particle) const;
-
- /**
- * Find an MCParticle from a Trajectory.
- * @param[in] trajectory The Trajectory.
- * @return The MCParticle or null if not found.
- */
- EVENT::MCParticle* findInitialParticle(Trajectory* trajectory) const;
-
- /**
- * Find a daughter MCParticle given a G4PrimaryParticle daughter.
- * @param[in] initial The parent particle.
- * @param[in] primary The G4PrimaryParticle daughter.
- * @return The primary daughter or null if does not exist.
- */
- EVENT::MCParticle* findDaughterParticle(EVENT::MCParticle* initial, G4PrimaryParticle* primary) const;
-
- /**
- * Find a final, output MCParticle given an initial, input MCParticle.
- * @param[in] initialParticle The intitial MCParticle.
- */
- IMPL::MCParticleImpl* findFinalParticle(EVENT::MCParticle* initialParticle) const;
-
- /**
- * Print out all the maps.
- */
- void printMaps();
-
- /**
- * Print the track to particle map.
- */
- void printTracksToParticles();
-
- /**
- * Print the particle to primary map.
- */
- void printParticlesToPrimaries();
-
- /**
- * Print the intial to final MCParticle map.
- */
- void printInitialParticleToFinal();
-
- /**
- * Clear the maps.
- */
- void clear();
-
-private:
-
- // trackID -> MCP
- TrackToMcpMap m_trackToMcp;
-
- // MCP -> G4PrimaryParticle
- McpToPrimaryMap m_McpToPrimary;
-
- // initial Mcp -> final Mcp
- InitMcpToMcpMap m_initMcpToMcp;
-
- // back pointer to LcioMcpManager
- LcioMcpManager* m_manager;
-};
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/LcioMcpMessenger.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/LcioMcpMessenger.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,59 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/LcioMcpMessenger.hh,v 1.12 2013-11-06 00:23:35 jeremy Exp $
-
-#ifndef SLIC_LCIOMCPMESSENGER_HH
-#define SLIC_LCIOMCPMESSENGER_HH 1
-
-// geant4
-#include "G4UIcommand.hh"
-#include "G4UIcmdWithADoubleAndUnit.hh"
-#include "G4UIcmdWithABool.hh"
-#include "G4UIdirectory.hh"
-#include "G4UImessenger.hh"
-
-namespace slic {
-
-/**
- * @class LcioMcpMessenger
- * @brief G4UImessenger for LcioMcpManager.
- */
-class LcioMcpMessenger: public G4UImessenger {
-
-public:
-
- /**
- * Class constructor.
- */
- LcioMcpMessenger();
-
- /**
- * Class destructor.
- */
- virtual ~LcioMcpMessenger();
-
-public:
-
- /**
- * Execute a macro command.
- * @param[in] command The command to execute.
- * @param[in] newValues The arguments.
- */
- virtual void SetNewValue(G4UIcommand* command, G4String newValues);
-
-private:
-
- /**
- * Define the commands handled by this messenger.
- */
- void defineCommands();
-
-private:
- G4UIdirectory* m_mcpDir;
-
- G4UIcmdWithADoubleAndUnit* m_minTrackingDistanceCmd;
- G4UIcmdWithABool* m_printFinalCollection;
- G4UIcmdWithABool* m_printInitialCollection;
- G4UIcmdWithABool* m_writeInitialCollection;
-};
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/LcioMcpPrinter.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/LcioMcpPrinter.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,165 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/LcioMcpPrinter.hh,v 1.21 2013-11-06 00:23:35 jeremy Exp $
-
-#ifndef SLIC_LCIOMCPPRINTER_HH
-#define SLIC_LCIOMCPPRINTER_HH 1
-
-// LCIO
-#include "EVENT/LCCollection.h"
-#include "EVENT/MCParticle.h"
-
-// SLIC
-#include "LcioMcpManager.hh"
-
-// STL
-#include <map>
-
-namespace slic {
-
-/**
- * @brief Prints MCParticle information for debugging.
- */
-class LcioMcpPrinter: public Module {
-public:
-
- /**
- * List of indices.
- */
- typedef std::vector<int> IndexVector;
-
-public:
-
- /**
- * Class constructor.
- */
- LcioMcpPrinter();
-
- /**
- * Class destructor.
- */
- virtual ~LcioMcpPrinter();
-
-public:
-
- /**
- * Print out an MCParticle collection.
- * @param[in] name The name of the collection.
- * @param[in] collection The LCCollection containing MCParticle objects.
- */
- void printParticleCollection(const std::string& name, EVENT::LCCollection* collection);
-
-private:
-
- /**
- * Print current MCParticle collection.
- */
- void printParticleCollection();
-
- /*
- * Print information for a single MCParticle.
- */
- void printParticle(EVENT::MCParticle* particle);
-
- /**
- * Make a string from a double array of length 3.
- * @param[in] array The double array.
- * @return The string.
- */
- static std::string makeDoubleArray3String(const double* array);
-
- /**
- * Make a string from a float array of length 3.
- * @param[in] array The float array.
- * @return The string.
- */
- static std::string makeFloatArray3String(const float* array);
-
- /**
- * Make a string from the status codes of an MCParticle.
- * @param[in] particle The MCParticle.
- * @return The list of status codes converted to a string.
- */
- static std::string makeSimStatusCodesString(EVENT::MCParticle* particle);
-
- /**
- * Make a string from a vector of indices.
- * @param[in] indices The IndexVector.
- * @return The list of indices converted to a string.
- */
- static std::string makeIndexVectorString(const IndexVector& indices);
-
- /**
- * Get the indices of an MCParticle's daughters.
- * @param[in] particle The MCParticle.
- * @return The vector of indices.
- */
- IndexVector findDaughterIndices(EVENT::MCParticle* particle);
-
- /**
- * Get the indices of an MCParticle's parents.
- * @param[in] particle The MCParticle.
- * @return The vector of indices.
- */
- IndexVector findParentIndices(EVENT::MCParticle* particle);
-
- /**
- * Find the index of a particle.
- * @param[in] particle The MCParticle.
- * @return The index of the MCParticle or -1 if not found.
- */
- int findParticleIndex(EVENT::MCParticle* particle);
-
- /**
- * Insert a field separator.
- */
- inline void sep();
-
- /**
- * Insert a new line.
- */
- inline void line();
-
- /**
- * Print the header.
- */
- inline void head();
-
- /**
- * Print the footer.
- */
- inline void foot();
-
- /**
- * Print the column labels.
- */
- void colLabels();
-
- /**
- * Print the status flag labels.
- */
- void flagLabels();
-
-public:
-
- static const int idx_width;
- static const int ptr_width;
- static const int pdg_width;
- static const int par_width;
- static const int dau_width;
- static const int mom_width;
- static const int vtx_width;
- static const int end_width;
- static const int time_width;
- static const int mass_width;
- static const int chrg_width;
- static const int e_width;
- static const int stat_width;
- static const int flags_width;
-
-private:
-
- EVENT::LCCollection* m_coll;
- std::string m_collName;
-};
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/LcioMcpStatusSetter.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/LcioMcpStatusSetter.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,78 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/LcioMcpStatusSetter.hh,v 1.7 2013-11-06 00:23:35 jeremy Exp $
-
-#ifndef SLIC_LCIOMCPSTATUSSETTER_HH
-#define SLIC_LCIOMCPSTATUSSETTER_HH 1
-
-// lcio
-#include "EVENT/MCParticle.h"
-#include "IMPL/MCParticleImpl.h"
-
-namespace slic {
-class Trajectory;
-
-/**
- * @class LcioMcpStatusSetter
- * @brief Static utility functions for setting the status of MCParticles
- * from associated Trajectory, initial MCParticle, and TrackInformation
- * objects.
- */
-class LcioMcpStatusSetter {
-
-private:
-
- /**
- * Class constructor, which is private because this class is all static utility methods.
- */
- LcioMcpStatusSetter() {
- }
-
-public:
-
- /**
- * Set the status codes on an MCParticle from a Trajectory.
- * @param[in] trajectory The particle Trajectory.
- * @param[in] particle The MCParticle on which to set status codes.
- * @return The same MCParticle with status codes set.
- */
- static IMPL::MCParticleImpl* setStatusCodesFromTrajectory(Trajectory* trajectory, IMPL::MCParticleImpl* particle);
-
- /**
- * Set the generator status on a final MCParticle from the associated initial.
- * @param[in] initialParticle The initial MCParticle.
- * @param[in,out] finalParticle The output MCParticle with generator status set.
- */
- static void setGeneratorStatus(EVENT::MCParticle* initialParticle, IMPL::MCParticleImpl* finalParticle);
-
-private:
-
- /**
- * Set the final flag values on an
- * @param[in] trajectory The Trajectory.
- * @param[in,out] particle The MCParticle on which statuses will be set.
- */
- static void setFinalStatusFlags(Trajectory* trajectory, IMPL::MCParticleImpl* particle);
-
- /**
- * Set the vertex is not endpoint of parent flag on an MCParticle from Trajectory data.
- * @param[in] trajectory The Trajectory.
- * @param[in,out] particle The MCParticle.
- */
- static void setVertexIsNotEndpointOfParentFlag(Trajectory* trajectory, IMPL::MCParticleImpl* particle);
-
- /**
- * Set the created in simulation flag on an MCParticle from Trajectory data.
- * @param[in] trajectory The Trajectory.
- * @param[in,out] particle The MCParticle.
- */
- static void setCreatedInSimulationFlag(Trajectory* trajectory, IMPL::MCParticleImpl* particle);
-
- /**
- * Set the backscatter flag on an MCParticle from Trajectory data.
- * @param[in] trajectory The Trajectory.
- * @param[in] particle The MCParticle.
- */
- static void setBackscatterFlag(Trajectory* trajectory, IMPL::MCParticleImpl* particle);
-};
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/LcioMcpUtil.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/LcioMcpUtil.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,54 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/LcioMcpUtil.hh,v 1.7 2013-11-06 00:23:35 jeremy Exp $
-
-#ifndef SLIC_LCIOMCPUTIL_HH
-#define SLIC_LCIOMCPUTIL_HH
-
-// lcio
-#include "IMPL/MCParticleImpl.h"
-
-// g4
-#include "globals.hh"
-
-namespace slic {
-
-/**
- * @class LcioMcpUtil
- * @brief Static MCParticle utility functions.
- */
-class LcioMcpUtil {
-
-private:
-
- /**
- * Class constructor, which is private because this class is all static utility functions.
- */
- LcioMcpUtil() {
- }
-
-public:
-
- /**
- * Compute the distance between the starting vertices of two MCParticle objects.
- * @param[in] particle1 The first MCParticle.
- * @param[in] particle2 The second MCParticle.
- * @return The distance in millimeters between the two particle's starting vertices.
- */
- static G4double computeParticleDistance(IMPL::MCParticleImpl* particle1, IMPL::MCParticleImpl* particle2);
-
- /**
- * Compute the distance between two points represented as double arrays.
- * @param[in] v1 The first vertex.
- * @param[in] v2 The second vertex.
- * @return The distance in millimeters between the two vertices.
- */
- static G4double computeDistance(const double *v1, const double *v2);
-
- /**
- * Check if an MCParticle is a primary, e.g. it has no parent particle.
- * @return True if MCParticle is primary; false if not.
- */
- static bool isPrimary(IMPL::MCParticleImpl* mcp);
-};
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/LcioMessenger.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/LcioMessenger.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -62,7 +62,8 @@
G4UIcommand* m_mcpFilenameCmd;
G4UIcmdWithAnInteger* m_setRunNumberCmd;
G4UIcommand* m_autonameCmd;
- G4UIcmdWithABool* m_dumpEventCmd;
+ G4UIcmdWithABool* m_printEventCmd;
+ G4UIcmdWithABool* _writeCompleteEventCmd;
};
}
projects/slic/trunk/include
--- projects/slic/trunk/include/LcioPrimaryGenerator.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/LcioPrimaryGenerator.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,76 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/LcioPrimaryGenerator.hh,v 1.18 2013-11-06 00:23:35 jeremy Exp $
-#ifndef SLIC_LCIOPRIMARYGENERATOR_HH
-#define SLIC_LCIOPRIMARYGENERATOR_HH
-
-// slic
-#include "Module.hh"
-
-// lcio
-#include "EVENT/LCCollection.h"
-#include "IMPL/MCParticleImpl.h"
-
-// geant4
-#include "G4VPrimaryGenerator.hh"
-#include "G4PrimaryParticle.hh"
-#include "G4PrimaryVertex.hh"
-#include "G4Event.hh"
-
-namespace slic {
-
-class LcioManager;
-class LcioMcpManager;
-
-/**
- @class LcioPrimaryGenerator
- @brief Generates G4PrimaryParticles from an LCIO MCParticle tree.
- */
-class LcioPrimaryGenerator: public Module {
-
-public:
- typedef std::map<EVENT::MCParticle*, G4PrimaryParticle*> ParticleToPrimaryMap;
-
-public:
-
- /**
- * Class constructor.
- * @param[in] manager The LcioManager.
- */
- LcioPrimaryGenerator(LcioManager* manager);
-
- /**
- * Class destructor.
- */
- virtual ~LcioPrimaryGenerator();
-
-public:
-
- /**
- * Create the G4Event structure from a collection of MCParticles.
- * @param[in] particles The list of MCParticles.
- * @param[in,out] anEvent The G4Event.
- */
- void generatePrimaryVertexFromMcpCollection(EVENT::LCCollection* particles, G4Event* anEvent);
-
-private:
-
- /**
- * Create the G4PrimaryParticle based on an input MCParticle.
- * @param[in] particle The input MCParticle.
- * @return The G4PrimaryParticle created from the MCParticle.
- */
- G4PrimaryParticle* createPrimaryParticle(IMPL::MCParticleImpl* particle);
-
- /**
- * Create a G4PrimaryVertex based on an itpu MCParticle.
- * @param[in] particle The MCParticle.
- * @return The G4PrimaryVertex created from the MCParticle.
- */
- G4PrimaryVertex* createPrimaryVertex(IMPL::MCParticleImpl* particle);
-
-private:
- LcioManager* m_mgr;
- LcioMcpManager* m_mcpManager;
-};
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/Module.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/Module.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -158,20 +158,6 @@
;
}
- //G4ClassificationOfNewTrack stackingClassifyNewTrack(const G4Track*) {
- // G4Exception("", "", FatalException, "This method should be overridden.");
- // Arbitrary return value. Will never reach here.
- // return fWaiting;
- //}
-
- //void stackingNewStage() {
- // ;
- //}
-
- //void stackingPrepareNewEvent() {
- // ;
- //}
-
public:
ModuleStatus m_status;
std::string m_name;
projects/slic/trunk/include
--- projects/slic/trunk/include/PackageInfo.hh.in 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/PackageInfo.hh.in 2014-07-15 23:59:54 UTC (rev 3187)
@@ -2,6 +2,8 @@
#ifndef SLIC_PACKAGE_INFO_HH
#define SLIC_PACKAGE_INFO_HH 1
+#include "lcdd/util/StringUtil.hh"
+
namespace slic {
/**
@@ -17,16 +19,27 @@
* Get the version string of the application.
* @return The version string.
*/
- static std::string& getVersionString() {
+ static std::string& getVersion() {
static std::string ver_str = "@SLIC_VERSION@";
return ver_str;
}
+
+ /**
+ * Get the version string with given seperator character and padding.
+ * @return The version string with seperator and padding.
+ */
+ static std::string getVersion(std::string& sep, int width = 2) {
+ std::string version = StringUtil::pad(std::string("@SLIC_VERSION_MAJOR@"), width);
+ version += sep + StringUtil::pad(std::string("@SLIC_VERSION_MINOR@"), width);
+ version += sep + StringUtil::pad(std::string("@SLIC_VERSION_PATCH@"), width);
+ return version;
+ }
/**
* Get the full name of the application.
* @return The full name of the application.
*/
- static std::string& getNameString() {
+ static std::string& getFullName() {
static std::string name_str = "Simulation for the Linear Collider";
return name_str;
}
@@ -35,8 +48,8 @@
* Get the abbreviated name of the application.
* @return The abbreviated name of the application.
*/
- static std::string& getAbbrevString() {
- static std::string abbrev_str = "slic";
+ static std::string& getShortName() {
+ static std::string abbrev_str = "SLIC";
return abbrev_str;
}
@@ -44,8 +57,8 @@
* Get a string fully describing the application.
* @return The description of the application.
*/
- static std::string getFullApplicationString() {
- static std::string fullVerStr = getNameString() + " " + getVersionString();
+ static std::string getApplication() {
+ static std::string fullVerStr = getShortName() + " " + getVersion();
return fullVerStr;
}
projects/slic/trunk/include
--- projects/slic/trunk/include/ParticleGunEventSource.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/ParticleGunEventSource.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -40,7 +40,7 @@
/**
* Dump information on the current event to the log.
*/
- void dumpCurrentEvent();
+ void printCurrentEvent();
/**
* Get the associated G4ParticleGun.
projects/slic/trunk/include
--- projects/slic/trunk/include/PrimaryGeneratorAction.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/PrimaryGeneratorAction.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -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/trunk/include
--- projects/slic/trunk/include/RunManager.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/RunManager.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -103,8 +103,6 @@
*/
void InitializeGeometry();
- //void DoEventLoop(G4int, const char*, G4int);
-
private:
/**
projects/slic/trunk/include
--- projects/slic/trunk/include/SlicApplication.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/SlicApplication.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -21,12 +21,20 @@
class SlicApplicationMessenger;
class RunManager;
-static const char* ErrorCodeStrings[] = { "okay", "unknown error", "caught unknown exception", "caught STD exception", "output file already exists",
- "failed to delete old output file", "event underflow", "invalid geometry setup", "event generator was never setup" };
+static const char* ErrorCodeStrings[] = {
+ "okay",
+ "unknown error",
+ "caught unknown exception",
+ "caught STD exception",
+ "output file already exists",
+ "failed to delete old output file",
+ "event underflow",
+ "invalid geometry setup",
+ "event generator was never setup" };
/**
* @class SlicApplication
- * @brief Singleton class for a Geant4 application.
+ * @brief Singleton class for a Geant4 application
* @note Instantiates manager and messenger classes, including
* user actions, G4RunManager, GeneratorManager, and PhysicsListManager.
*/
@@ -38,7 +46,8 @@
* Indicates the run mode: batch or interactive.
*/
enum ERunMode {
- eBatch = 1, eInteractive = 2
+ eBatch = 1,
+ eInteractive = 2
};
/**
@@ -169,12 +178,6 @@
void setReturnCode(int rc);
/**
- * Get the Geant4 version string.
- * @return The Geant4 version string.
- */
- std::string getGeant4VersionString();
-
- /**
* Get a string representation of an error code.
* @param[in] errorCode The errorCode.
* @return The string representation of the error code.
@@ -230,9 +233,6 @@
// return code
int m_returnCode;
- // hack to make aborting G4 work (doesn't seem to function properly)
- bool m_setRunAbort;
-
// Has the initialize() function been called?
bool m_isInitialized;
projects/slic/trunk/include
--- projects/slic/trunk/include/StdHepEventSource.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/StdHepEventSource.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -3,10 +3,9 @@
#ifndef SLIC_STDHEPEVENTSOURCE_HH
#define SLIC_STDHEPEVENTSOURCE_HH 1
-// slic
-#include "StdHepLoader.hh"
-#include "StdHepToLcioConvertor.hh"
+// SLIC
#include "EventSourceWithInputFile.hh"
+#include "StdHepGenerator.hh"
namespace slic {
/**
@@ -50,7 +49,7 @@
/**
* Print out the current event.
*/
- void dumpCurrentEvent();
+ void printCurrentEvent();
/**
* Begin of run action.
@@ -64,15 +63,8 @@
*/
void beginEvent(const G4Event* anEvent);
- /**
- * Get the lStdHep reader associated with this source.
- * @return The StdHep reader.
- */
- lStdHep* getStdHepReader();
-
private:
- StdHepLoader* m_loader;
- StdHepToLcioConvertor* m_convertor;
+ StdHepGenerator* m_eventGenerator;
};
}
projects/slic/trunk/include
--- projects/slic/trunk/include/StdHepLoader.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/StdHepLoader.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,98 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/StdHepLoader.hh,v 1.18 2013-11-06 00:23:35 jeremy Exp $
-
-#ifndef slic_StdhepLoader_hh
-#define slic_StdhepLoader_hh 1
-
-// slic
-#include "Module.hh"
-
-// lcio
-#include "UTIL/lStdHep.hh"
-
-// std
-#include <string>
-#include <fstream>
-
-class G4Event;
-class G4PrimaryVertex;
-class G4PrimaryParticle;
-
-using UTIL::lStdHep;
-
-namespace slic {
-enum EErrorCode {
- eEOF = 106
-};
-
-/**
- @class StdHepLoader
- @brief Simple wrapper to lStdHep interface.
- */
-class StdHepLoader: public Module {
-
-public:
-
- /**
- * Class constructor.
- */
- StdHepLoader();
-
- /**
- * Class destructor.
- */
- virtual ~StdHepLoader();
-
-public:
-
- /**
- * Get the file name.
- * @return The file name.
- */
- const std::string& getFilename() const {
- return m_filename;
- }
-
- /**
- * Get the backing lStdHep reader.
- * @return The lStdHep reader object.
- */
- lStdHep* getStdHepReader() {
- return m_reader;
- }
-
- /**
- * Return true if end of input events.
- * @return True if end of input events; false if not.
- */
- bool isEndOfInput() {
- return m_isEndOfInput;
- }
-
-public:
-
- /**
- * Open a StdHep file for reading.
- * @param[in] filename The full path to the file.
- */
- void openStdHepFile(const std::string& filename);
-
- /**
- * Open a StdHep file for reading.
- * @param[in] filename The full path to the file.
- */
- void openStdHepFile(const char* filename);
-
- void closeStdHepFile();
-
- void readNextEvent();
- void dumpCurrentEvent();
-
-private:
-
- std::string m_filename;
- lStdHep* m_reader;
- bool m_isEndOfInput;
-};
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/StdHepToLcioConvertor.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/StdHepToLcioConvertor.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,195 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/StdHepToLcioConvertor.hh,v 1.23 2013-11-06 00:23:35 jeremy Exp $
-
-#ifndef SLIC_STDHEPTOLCIOCONVERTOR_HH
-#define SLIC_STDHEPTOLCIOCONVERTOR_HH 1
-
-// lcdd
-#include "Module.hh"
-
-// lcio
-#include "IMPL/MCParticleImpl.h"
-#include "EVENT/LCCollection.h"
-#include "UTIL/lStdHep.hh"
-
-// std
-#include <vector>
-#include <cassert>
-
-using UTIL::lStdHep;
-
-namespace slic {
-
-/**
- * @class StdHepToLcioConvertor
- * @brief Converts an input StdHep event to an LCIO MCParticle collection.
- */
-
-/**
- @verbatim
-
- The index "pointers" daughter1, daughter2, mother1 and mother2 are handled as follows,
- where idx1 is daughter1/mother1, and idx2 is daughter2/mother2.
-
- idx1 cmp idx2 action
- _______________________
-
- -1 -1 do nothing
- # -1 use idx1
- -1 # use idx2
- # == # use idx1
- # < # range = idx1 to idx2
- # > # use idx1 and idx2 as discrete indices
-
- @endverbatim
- */
-class StdHepToLcioConvertor: public Module {
-
-public:
-
- /**
- * Class constructor.
- * @param[in] reader The lStdHep to use.
- */
- StdHepToLcioConvertor(lStdHep* reader);
-
- /**
- * Class constructor.
- */
- StdHepToLcioConvertor();
-
- /**
- * Class destructor.
- */
- virtual ~StdHepToLcioConvertor();
-
-public:
-
- /**
- * Set the lStdHep reader to be used by the convertor.
- * @param[in] rdr The lStdHep reader to use.
- */
- void setStdHepReader(lStdHep* rdr) {
- assert(rdr);
- m_reader = rdr;
- }
-
- /**
- * Create the next MCP collection using current event
- * @note This function depends upon an external read operation to setup the current event.
- */
- void fillParticleCollection();
-
- /**
- * Reset the current MCParticle collection.
- * @param[in] collection The new MCParticle collection.
- */
- void setParticleCollection(EVENT::LCCollection* collection);
-
- /**
- * Get the current MCParticle collection.
- * @return The current MCParticle collection.
- */
- EVENT::LCCollection* getParticleCollection();
-
- /**
- * Create a single MCParticle from a StdHep record.
- * @param[in] ihep The index of the StdHep record.
- * @note Does not setup mother or daughter links.
- *
- */
- IMPL::MCParticleImpl* createParticle(int ihep);
-
- /**
- * Setup the parent pointers of this MCParticle.
- * @param[in] ihep The index of the StdHep record.
- * @param[out] particle The MCParticle which will have its parent pointers setup.
- */
- void setupParents(int ihep, IMPL::MCParticleImpl* particle);
-
- /**
- * Assign parent MCParticle to the daugher from the index.
- * @param[in] parentIndex The MCParticle index of the parent.
- * @param[out] daughter The daughter particle.
- */
- void addParent(int parentIndex, IMPL::MCParticleImpl* daughter);
-
- /**
- * Assign the StdHep record at the index as a daughter track of the MCParticle.
- * @param[in] ihep The StdHep record index.
- * @param[out] particle The parent particle.
- */
- void setupDaughters(int ihep, IMPL::MCParticleImpl* particle);
-
- /**
- * Setup the input MCParticle as a parent of this track.
- * @param[in] daughterIndex The index of the daughter.
- * @param[in] parent The parent MCParticle.
- */
- void addDaughter(int daughterIndex, IMPL::MCParticleImpl* parent);
-
- /**
- * Check if an MCParticle has another particle as its parent.
- * @param[in] daughter The proposed daughter MCParticle.
- * @param[in] parent The parent MCParticle.
- * @return True if daughter particle has this parent; false if not.
- */
- bool hasParent(IMPL::MCParticleImpl* daughter, IMPL::MCParticleImpl* parent);
-
- /**
- * Given two ihep indices, return a vector of all the included indices.
- * @param[in] startIndex The starting index.
- * @param[in] endIndex The ending index.
- * @param[out] vec The vector of indices.
- */
- std::vector<int> makeIndexVec(int startIndex, int endIndex, std::vector<int>& vec);
-
- /**
- * Perform various sanity checks on the track's mother and daughter index values.
- * @param[in] The StdHep record's index.
- */
- void checkParentage(int ihep);
-
- /**
- * Print the index.
- * @param[in] ihep The index.
- */
- void printIndex(int ihep);
-
- /**
- * Print the mothers of this track.
- * @param[in] ihep The index.
- */
- void printMothers(int ihep);
-
- /**
- * Print the daughters of this track.
- * @param[in] ihep The index.
- */
- void printDaughters(int ihep);
-
- /**
- * Print the track.
- * @param[in] The index.
- */
- void printTrack(int ihep);
-
- /**
- * Return true if index is basically valid, in that it is greater than zero
- * and not greater than the total number of tracks in the event.
- * @return True if daughter index is valid; false if not.
- */
- bool isValidIndex(int ihep) {
- return (ihep > -1 && ihep < m_reader->nTracks());
- }
-
-private:
-
- /* Pointer to the current MCParticle collection that is being filled. */
- EVENT::LCCollection* m_currentMcpColl;
-
- /* Pointer to the lStdHep reader being used as a data source. */
- lStdHep* m_reader;
-};
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/SteppingAction.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/SteppingAction.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -40,6 +40,13 @@
*/
void UserSteppingAction(const G4Step* step);
+ /**
+ * Return true if this track appears to be backscattering from a non-tracking region
+ * into the tracking region.
+ * @return True if track is backscattering; false if not.
+ */
+ bool isBackScattering(const G4Step* step);
+
private:
/**
projects/slic/trunk/include
--- projects/slic/trunk/include/TrackingAction.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/TrackingAction.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -12,7 +12,7 @@
/**
* @class TrackingAction
- * @brief Implementation of G4UserTrackingAction.
+ * @brief Implementation of G4UserTrackingAction
*/
class TrackingAction: public G4UserTrackingAction {
@@ -21,12 +21,14 @@
/**
* Class constructor.
*/
- TrackingAction();
+ TrackingAction() {
+ }
/**
* Class destructor.
*/
- ~TrackingAction();
+ ~TrackingAction() {
+ }
public:
@@ -34,17 +36,13 @@
* Pre tracking action.
* @param[in] aTrack The G4Track.
*/
- void PreUserTrackingAction(const G4Track* aTrack);
+ void PreUserTrackingAction(const G4Track*);
/**
* Post tracking action.
* @param[in] aTrack The G4Track.
*/
- void PostUserTrackingAction(const G4Track* aTrack);
-
-private:
-
- TrajectoryManager* m_trjMgr;
+ void PostUserTrackingAction(const G4Track*);
};
}
projects/slic/trunk/include
--- projects/slic/trunk/include/Trajectory.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/Trajectory.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,337 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/Trajectory.hh,v 1.19 2013-11-06 00:23:36 jeremy Exp $
-
-#ifndef SLIC_TRAJECTORY_HH
-#define SLIC_TRAJECTORY_HH 1
-
-// geant4
-#include "G4VTrajectory.hh"
-#include "G4Allocator.hh"
-#include "G4Track.hh"
-#include "G4VTrajectoryPoint.hh"
-#include "globals.hh"
-
-namespace slic {
-
-typedef std::vector<G4VTrajectoryPoint*> TrajectoryPointContainer;
-
-/**
- * @class Trajectory
- * @brief Implementation of G4VTrajectory.
- */
-class Trajectory: public G4VTrajectory {
-public:
-
- /**
- * Final status flags.
- */
- enum EFinalStatus {
- eUnset = 0, eDecayedInTracker = 1, eLeftDetector = 2, eStopped = 3, eDecayedInCalorimeter = 4
- };
-
-public:
-
- /**
- * Class constructor.
- */
- Trajectory();
-
- /**
- * Class constructor.
- * @param[in] aTrack An associated G4Track.
- */
- Trajectory(const G4Track *aTrack);
-
- /**
- * Copy constructor.
- * @param[in] trajectory The Trajectory to copy.
- */
- Trajectory(Trajectory& trajectory);
-
- /**
- * Class destructor.
- */
- virtual ~Trajectory();
-
-public:
-
- /**
- * Draw a Trajectory.
- * @param[in] i_mode The draw mode.
- */
- //virtual void DrawTrajectory(G4int i_mode) const {
- //
- virtual void DrawTrajectory() const {
- G4VTrajectory::DrawTrajectory();
- }
-
- /**
- * Convert final status enum to string (static version).
- * @param[in] fs The EFinalStatus.
- * @return The status converted to a string.
- */
- static const std::string& getFinalStatusString(EFinalStatus fs);
-
- /**
- * Convert the current final status enum to string.
- * @return The current final status converted to a string.
- */
- inline const std::string& getFinalStatusString() const {
- return Trajectory::getFinalStatusString(m_finalStatus);
- }
-
- /**
- * Specialized implementation of the new operator.
- * @param[in] size The size to allocate.
- */
- inline void* operator new(size_t size);
-
- /**
- * Specialized implementation of the delete operator.
- * @param[in] pointer Pointer to the object to delete.
- */
- inline void operator delete(void* pointer);
-
- /**
- * Equals operator.
- * @param[in] right The comparison object.
- */
- inline int operator ==(const Trajectory& right) const {
- return (this == &right);
- }
-
- /**
- * Get the track ID of this Trajectory.
- * @return The track ID.
- */
- inline G4int GetTrackID() const {
- return m_trackID;
- }
-
- /**
- * Get parent track ID.
- * @return The parent track ID.
- */
- inline G4int GetParentID() const {
- return m_parentID;
- }
-
- /**
- * Get the name of the particle associated with the Trajectory.
- * @return The name of the particle.
- */
- inline G4String GetParticleName() const {
- return m_particleName;
- }
-
- /**
- * Get the charge.
- * @return The charge.
- */
- inline G4double GetCharge() const {
- return m_PDGCharge;
- }
-
- /**
- * Get the PDG encoding of this Trajectory.
- * @return The PDG encoding.
- */
- inline G4int GetPDGEncoding() const {
- return m_PDGEncoding;
- }
-
- /**
- * Get the initial momentum.
- * @return The initial momentum.
- */
- inline G4ThreeVector GetInitialMomentum() const {
- return m_initialMomentum;
- }
-
- /**
- * Get the particle definition.
- * @return The particle definition.
- */
- inline const G4ParticleDefinition* GetParticleDefinition() {
- return m_particleDefinition;
- }
-
- /**
- * Get the creation time of the particle, measured from the start of the event.
- * @return The creation time.
- */
- inline G4double GetGlobalTime() {
- return m_globalTime;
- }
-
- /**
- * Return true if particle has daughters produced at its endpoint.
- * @return True if particle has endpoint daughters; false if not.
- */
- inline G4bool getHasEndpointDaughters() const {
- return m_hasEndpointDaughters;
- }
-
- /**
- * Set that this Trajectory has endpoint daughters.
- * @param[in] d Set to true if has endpoint daughters; false if not.
- */
- inline void setHasEndpointDaughters(G4bool d) {
- m_hasEndpointDaughters = d;
- }
-
- /**
- * Return true if this Trajectory is backscattered from a non-tracking region
- * back into the tracking region.
- * @return True if Trajectory is backscattered; false if not.
- */
- inline G4bool getBackscatter() const {
- return m_backscatter;
- }
-
- /**
- * Set backscatter.
- * @param[in] b Set to true if backscatter; false if not.
- */
- inline void setBackscatter(G4bool b) {
- m_backscatter = b;
- }
-
- /**
- * Get vertex is not endpoint of parent flag.
- * @return The vertix is not endpoint of parent flag.
- */
- inline G4bool getVertexIsNotEndpointOfParent() const {
- return m_vertexIsNotEndpointOfParent;
- }
-
- /**
- * Set the vertex is not endpoint of parent flag.
- * @param[in] v The vertex is not endpoint of parent flag value.
- */
- inline void setVertexIsNotEndpointOfParent(G4bool v = true) {
- m_vertexIsNotEndpointOfParent = true;
- }
-
- /**
- * Set the final status.
- * @param[in] fstatus The final status.
- */
- inline void setFinalStatus(EFinalStatus fstatus) {
- m_finalStatus = fstatus;
- }
-
- /**
- * Get the final status.
- * @return The final status.
- */
- inline EFinalStatus getFinalStatus() const {
- return m_finalStatus;
- }
-
- /**
- * Set the created in simulation flag.
- * @param[in] c The created in simulation flag.
- */
- inline void setCreatedInSimulation(G4bool c) {
- m_createdInSimulation = c;
- }
-
- /**
- * Get the created in simulation flag.
- * @return The created in simulation flag.
- */
- inline G4bool getCreatedInSimulation() const {
- return m_createdInSimulation;
- }
-
- /**
- * Setup this Trajectory from a G4Track.
- * @param[in] aTrack The G4Track.
- */
- void setupManuallyFromTrack(const G4Track* aTrack);
-
- /**
- * Set the endpoint energy.
- * @param[in] e The endpoint energy.
- */
- inline void setEndpointEnergy(G4double e) {
- m_endpointEnergy = e;
- }
-
- /**
- * Get the endpoint energy.
- * @return The endpoint energy in GeV.
- */
- G4double getEndPointEnergy() {
- return m_endpointEnergy;
- }
-
-public:
-
- /**
- * Append the information from a G4Step.
- * @param[in] aStep The G4Step.
- */
- void AppendStep(const G4Step *aStep);
-
- /**
- * Merge a trajectory into this one.
- * @param[in] secondTrajectory The Trajectory to be merged.
- */
- void MergeTrajectory(G4VTrajectory *secondTrajectory);
-
- /**
- * Get the number of point entries.
- * @return The number of point entries.
- */
- virtual G4int GetPointEntries() const {
- return m_posRecord->size();
- }
-
- /**
- * Get a G4VTrajectoryPoint by index.
- * @param[in] i The index.
- * @return The trajectory point at index.
- */
- virtual G4VTrajectoryPoint* GetPoint(G4int i) const {
- return (*m_posRecord)[i];
- }
-
-private:
-
- // track info
- G4ParticleDefinition* m_particleDefinition;
- G4String m_particleName;
- G4double m_PDGCharge;
- G4int m_PDGEncoding;
- G4int m_trackID;
- G4int m_parentID;
- TrajectoryPointContainer* m_posRecord;
- G4ThreeVector m_initialMomentum;
- G4double m_globalTime;
- G4double m_endpointEnergy;
-
- // non-exclusive sim statuses
- G4bool m_createdInSimulation;
- G4bool m_hasEndpointDaughters;
- G4bool m_backscatter;
- G4bool m_vertexIsNotEndpointOfParent;
-
- // mutually exclusive final states: see enums
- EFinalStatus m_finalStatus;
-};
-
-extern G4Allocator<Trajectory> TrajectoryAllocator;
-
-inline void* Trajectory::operator new(size_t) {
- void *aTrajectory;
- aTrajectory = (void*) TrajectoryAllocator.MallocSingle();
- return aTrajectory;
-}
-
-inline void Trajectory::operator delete(void* aTrajectory) {
- TrajectoryAllocator.FreeSingle((Trajectory*) aTrajectory);
-}
-}
-
-#endif
projects/slic/trunk/include
--- projects/slic/trunk/include/TrajectoryManager.hh 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/include/TrajectoryManager.hh 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,405 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/include/TrajectoryManager.hh,v 1.36 2013-11-06 00:23:36 jeremy Exp $
-
-#ifndef SLIC_TRAJECTORYMANAGER_HH
-#define SLIC_TRAJECTORYMANAGER_HH 1
-
-// LCDD
-#include "lcdd/geant4/UserRegionInformation.hh"
-#include "lcdd/hits/TrackInformation.hh"
-
-// SLIC
-#include "Module.hh"
-#include "Singleton.hh"
-#include "Trajectory.hh"
-
-// Geant4
-#include "G4Track.hh"
-#include "G4TrackingManager.hh"
-
-// STL
-#include <map>
-
-class G4Event;
-
-namespace slic {
-
-/**
- * @class TrajectoryManager
- *
- * @brief Provides runtime management of Trajectory and TrackInformation objects,
- * mainly to update them so that the LcioMcpManager can apply the correct status
- * settings to output MCParticles.
- *
- * @todo Add current track ID and G4Track instance variables to simplify
- * internal handling and function signatures.
- */
-
-class TrajectoryManager: public Module, public Singleton<TrajectoryManager> {
-
-public:
-
- /**
- * Map of track ID to Trajectory.
- */
- typedef std::map<G4int, Trajectory*> TrackIDToTrajectoryMap;
-
-public:
-
- /**
- * Class destructor.
- */
- virtual ~TrajectoryManager();
-
- /**
- * Class constructor.
- */
- TrajectoryManager();
-
-public:
-
- /**
- * Begin of event action.
- * Clears the TrackIDToTrajectoryMap for new data.
- * @param[in] event The G4Event that is beginning.
- */
- void beginEvent(const G4Event* event) {
- clearTrackIDToTrajectoryMap();
- }
-
- /**
- * Pre-tracking action.
- * @param[in] track The new G4Track.
- */
- void preTracking(const G4Track* track);
-
- /**
- * Post-tracking action.
- * @param[in] track The G4Track that is about to be deleted.
- */
- void postTracking(const G4Track* track);
-
- /**
- * The stepping action.
- * @param[in] step The current G4Step.
- */
- void stepping(const G4Step* step);
-
- /**
- * Find a Trajectory by trackID.
- * @param[in] trackID The track ID.
- * @return The matching Trajectory or null if not found.
- */
- Trajectory* findTrajectory(G4int trackID);
-
- /**
- * Find a Trajectory from a G4Track.
- * @param[in] aTrack The G4Track.
- * @return The matching Trajectory or null if not found.
- */
- Trajectory* findTrajectory(const G4Track* aTrack) {
- return findTrajectory(aTrack->GetTrackID());
- }
-
- /**
- * Find a Trajectory from a G4Track.
- * @param[in] aTrack The G4Track.
- * @return The matching Trajectory or null if not found.
- */
- Trajectory* findTrajectory(G4Track* aTrack) {
- return findTrajectory((const_cast<const G4Track*>(aTrack))->GetTrackID());
- }
-
- /**
- * Print out the TrackIDToTrajectoryMap.
- */
- void printTrackIDToTrajectoryMap() {
- log() << LOG::debug << "TrackIDToTrajectoryMap" << LOG::done;
- log() << LOG::debug << "trackID | trjPtr" << LOG::done;
- for (TrackIDToTrajectoryMap::iterator iter = m_trackIDToTrajectory.begin(); iter != m_trackIDToTrajectory.end(); iter++) {
- log() << LOG::debug << iter->first << " " << iter->second << LOG::done;
- }
- }
-
- /**
- * Set the G4TrackingManager.
- * @param[in] trackManager The G4TrackingManager to be used.
- */
- inline void setTrackingManager(const G4TrackingManager* trackManager) {
- m_trackingManager = const_cast<G4TrackingManager*>(trackManager);
- }
-
-private:
-
- /**
- * Handle a primary track, e.g. one without a parent.
- * This function will set its TrackInformation and create a Trajectory for it.
- * @param[in] The G4Track which is a primary.
- */
- void handlePrimaryTrack(const G4Track* aTrack);
-
- /**
- * Setup the TrackInformation for a primary track.
- * @param[in] aTrack The G4Track which is a primary.
- */
- void setPrimaryTrackInformation(const G4Track* aTrack);
-
- /**
- * Handle a secondary track, which is the daughter of another.
- * This method will set the status flags on this track's Trajectory,
- * or create a TrackInformation object and Trajectory for it,
- * if needed.
- * @param[in] aTrack The G4Track which is a secondary.
- */
- void handleSecondaryTrack(const G4Track* aTrack);
-
- /**
- * Create a Trajectory for a secondary G4Track. Called from handleSecondaryTrack().
- * @param[in] aTrack The secondary G4Track.
- */
- Trajectory* createSecondaryTrajectory(const G4Track* aTrack);
-
- /**
- * Setup the TrackInformation for a secondary.
- * @param[in] aTrack The secondary G4Track.
- */
- void setupSecondaryTrackInformation(const G4Track* aTrack);
-
- /**
- * Copy the G4Track's TrackInformation to its secodaries.
- * @param[in] aTrack The parent G4Track.
- */
- void copySecondaryTrackInformationFromParent(const G4Track* aTrack);
-
- /**
- * Check if this track requires manual trajectory setup, e.g. if it
- * has a SimTrackerHit assigned to it but no Trajectory is being stored.
- * @param[in] aTrack The G4Track.
- * @return True if G4Track needs a Trajectory to be setup; false if not.
- */
- bool needsManualTrajectorySetup(const G4Track* aTrack);
-
- /**
- * Set the final status flag of a G4Track.
- * @param[in] aTrack The G4Track.
- */
- void setTrajectoryFinalStatus(const G4Track* aTrack);
-
- /**
- * Get the number of secondaries at the endpoint of the track.
- * @param[in] aTrack The G4Track.
- * @return The number of secondaries produced at the end of the track.
- */
- G4int getNumberOfSecondariesAtEnd(const G4Track* aTrack);
-
- /**
- * Return true if this track is primary, i.e. it has no parent.
- * @return True if this track is primary; false if it is a secondary.
- */
- inline bool isPrimaryTrack(const G4Track* aTrack) const {
- return (aTrack->GetParentID() == 0);
- }
-
- /**
- * Return true if energy is less than region's threshold to create a trajectory or if this
- * flag has been set on the TrackInformation.
- * @return True if track is below threshold; false if not.
- */
- inline bool isBelowThreshold(const G4Track* aTrack) const {
- return ((aTrack->GetKineticEnergy() < UserRegionInformation::getRegionInformation(aTrack)->getThreshold())
- || (TrackInformation::getTrackInformation(aTrack)->getBelowThreshold() == true));
- }
-
- /**
- * Check if the particle is backscattered back into a tracking region.
- * This checks that the original tracking status is "not in a tracking region",
- * which works ok because this function is only relevent when the track is
- * currently within the tracking region.
- * @return true if track is backscattered; false if not.
- */
- inline bool isBackscatter(const G4Track* aTrack) const {
- return TrackInformation::getTrackInformation(aTrack)->getOriginalTrackingStatus() != TrackInformation::eInTrackingRegion;
- }
-
- /**
- * Return true if the G4Track associated with the G4Step is in the "alive" state.
- * @return true if the G4Track from the G4Step is alive; false if not.
- */
- bool isAlive(const G4Step* aStep) {
- return aStep->GetTrack()->GetTrackStatus() == fAlive;
- }
-
- /**
- * Get the current Trajectory from the G4TrackingManager.
- * @return The Trajectory from the G4TrackingManager.
- */
- Trajectory* getCurrentTrajectory() {
- return static_cast<Trajectory*>(m_trackingManager->GimmeTrajectory());
- }
-
- /**
- * Fill the track ID to Trajectory map from the G4Event.
- * @param[in] event The G4Event.
- * @param[in] clearIt Set to true to clear the map before filling it.
- */
- void fillTrajectoryMap(const G4Event* event, bool clearIt = true);
-
- /**
- * Add a Trajectory to map by its track ID.
- * @param[in] trajectory The Trajectory to add.
- */
- void addTrackIDToTrajectoryLink(Trajectory* trajectory);
-
- /**
- * Clear the TrackIDToTrajectoryMap.
- */
- void clearTrackIDToTrajectoryMap() {
- m_trackIDToTrajectory.clear();
- }
-
- /**
- * Set the tracking flag from the UserRegionInformation, i.e. in or outside tracking region.
- * @param[out] trackInformation The TrackInformation to update.
- * @param[in] regionInformation The region information with the status.
- */
- void setTrackingFlagFromRegionInformation(TrackInformation* trackInformation, UserRegionInformation* regionInformation);
-
- /**
- * Using the parent's endpoint and the secondaries starting position to determine
- * the value of the vertexIsNotEndpointOfParent flag.
- * @param[in] endpoint The parent's endpoint.
- * @param[in] secondary The secondary G4Track.
- * @return true if vertex is not endpoint of parent; false if not.
- */
- bool vertexIsNotEndpointOfParent(const G4ThreeVector& endpoint, G4Track* secondary) {
- //assert(secondary);
-
- bool r = false;
-
- if (!endpoint.isNear(secondary->GetPosition())) {
- r = true;
- }
-
- return r;
- }
-
- /**
- * Get the Trajectory's point position by index.
- * @param[in] trajectory The Trajectory.
- * @param[in] positionIndex The index into the vector of point positions.
- * @return The Trajectory's point position.
- */
- static G4ThreeVector getTrajectoryPointPosition(Trajectory* trajectory, G4int positionIndex = 0) {
- // check null trajectory
- if (!trajectory) {
- G4Exception("TrajectoryManager::getTrajectoryPointPosition()", "", FatalException, "Trajectory is null.");
- }
-
- // get number of trajectories
- G4int numPnts = trajectory->GetPointEntries();
-
- // check idx argument against number of trajectories
- if ((positionIndex > (numPnts - 1)) || positionIndex < 0) {
- G4Exception("TrajectoryManager::getTrajectoryPointPosition()", "", FatalException, "The trajectory position index is invalid.");
- }
-
- // return point at idx
- return trajectory->GetPoint(positionIndex)->GetPosition();
- }
-
- /**
- * Get a Trajectory's vertex point, which is the first entry in the point position vector.
- * @param[in] The Trajectory.
- * @return The point position at the vertex.
- */
- inline static G4ThreeVector getTrajectoryVertex(Trajectory* trajectory) {
- return getTrajectoryPointPosition(trajectory, 0);
- }
-
- /**
- * Get a trajectory's endpoint, which is the last entry in the point position vector.
- * @param[in] trajectory The Trajectory.
- * @return The Trajectory's endpoint position.
- */
- inline static G4ThreeVector getTrajectoryEndpoint(Trajectory* trajectory) {
- return getTrajectoryPointPosition(trajectory, trajectory->GetPointEntries() - 1);
- }
-
- /**
- * Create a new trajectory from a G4Track.
- * @param[in] aTrack The source G4Track.
- * @return The output Trajectory from the track data.
- */
- Trajectory* createTrajectoryFromTrack(const G4Track* aTrack);
-
- /**
- * Get the current Trajectory being processed.
- * @return Pointer to the current trajectory being processed.
- */
- inline Trajectory* getCurrentTrajectory() const {
- return m_currentTrajectory;
- }
-
- /**
- * Get the current G4Track being processed.
- * @return The current G4Track.
- */
- inline const G4Track* getCurrentTrack() const {
- return m_currentTrack;
- }
-
- /**
- * Get the current track ID.
- * @return The current G4Track's ID.
- */
- inline G4int getCurrentTrackID() const {
- return m_currentTrackID;
- }
-
- /**
- * Print the status flags of a G4Track.
- * @param[in] track The G4Track.
- */
- void printStatusFlags(const G4Track* track);
-
- /**
- * Print the status flags from the TrackInformation.
- * @param[in] trackInformation The TrackInformation.
- */
- void printStatusFlags(TrackInformation* trackInformation);
-
- /**
- * Print the data of the G4Track.
- * @param[in] track The G4Track.
- * @param[in] isSecondary Set to true if track is a secondary; false if not.
- */
- void printTrack(const G4Track* track, bool isSecondary = false);
-
- /**
- * Print out the Trajectory.
- * @param[in] trajectory The Trajectory.
- */
- void printTrajectory(const Trajectory* trajectory);
-
- /**
- * Print out the secondaries of the current track.
- */
- void printSecondaries();
-
- /**
- * Print out physical volume information from the G4StepPoint.
- * @param[in] stepPoint The G4StepPoint.
- */
- void printPhysVolDbg(G4StepPoint* stepPoint);
-
-
-private:
-
- G4TrackingManager* m_trackingManager;
- TrackIDToTrajectoryMap m_trackIDToTrajectory;
- Trajectory* m_currentTrajectory;
- const G4Track* m_currentTrack;
- G4int m_currentTrackID;
-};
-}
-
-#endif
projects/slic/trunk/macros
--- projects/slic/trunk/macros/gps.mac 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/macros/gps.mac 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,3 +1,4 @@
+/generator/select gps
/gps/verbose 2
/gps/particle mu-
/gps/energy 10 GeV
projects/slic/trunk/src
--- projects/slic/trunk/src/EventAction.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/EventAction.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -5,12 +5,13 @@
#include "lcdd/util/StringUtil.hh"
// SLIC
-#include "LcioManager.hh"
-#include "TrajectoryManager.hh"
#include "EventSourceManager.hh"
-#include "SlicApplication.hh"
#include "EventMessenger.hh"
-#include "EventDebugger.hh"
+#include "LcioManager.hh"
+#include "MCParticleManager.hh"
+#include "RunManager.hh"
+//#include "SlicApplication.hh"
+#include "TrackManager.hh"
// Geant4
#include "G4RunManager.hh"
@@ -28,48 +29,35 @@
}
void EventAction::BeginOfEventAction(const G4Event *anEvent) {
- // check aborted
- if (!anEvent->IsAborted()) {
+}
- // debugger
- //EventDebugger::instance()->beginEvent( anEvent );
+void EventAction::EndOfEventAction(const G4Event *anEvent) {
- // trj mgr
- TrajectoryManager::instance()->beginEvent(anEvent);
+ /* Only execute these hooks if run was not aborted. */
+ if (!RunManager::instance()->isRunAborted()) {
- /* MCP begin event */
- LcioMcpManager::instance()->beginEvent(anEvent);
- }
-#ifdef SLIC_LOG
- else {
- log().warning("aborted EventAction::BeginOfEventAction");
- }
-#endif
-}
+ /* Create an empty LCEvent. */
+ LcioManager::instance()->createLCEvent();
-void EventAction::EndOfEventAction(const G4Event *anEvent) {
- // check aborted
- if (!SlicApplication::instance()->isAborting()) {
- // LcioManager's action
+ /* Save track information to MCParticle collection in LCIO output event. */
+ TrackManager::instance()->saveTrackSummaries(anEvent, LcioManager::instance()->getCurrentLCEvent());
+
+ /* End of event processing for the current event generator. */
+ EventSourceManager::instance()->endEvent(anEvent);
+
+ /* Execute LcioManager's end of event action. */
LcioManager::instance()->endEvent(anEvent);
}
- // event source (generator) action
- EventSourceManager::instance()->endEvent(anEvent);
-
- // event timer
+ /* Stop the event timer. */
stopEventTimer();
- // debugger
- //EventDebugger::instance()->endEvent(anEvent);
-
- // end event mesg
+ /* Print the end event message. */
printEndEventMessage(anEvent);
}
void EventAction::printEndEventMessage(const G4Event *anEvent) {
- log() << LOG::okay << ">>>> EndEvent <" + StringUtil::toString(anEvent->GetEventID()) + ">" << LOG::endl
- << LOG::done;
+ log() << LOG::okay << ">>>> EndEvent <" + StringUtil::toString(anEvent->GetEventID()) + ">" << LOG::endl << LOG::done;
}
EventAction* EventAction::getEventAction() {
projects/slic/trunk/src
--- projects/slic/trunk/src/EventDebugger.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/EventDebugger.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,113 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/EventDebugger.cc,v 1.7 2012-11-27 19:32:19 jeremy Exp $
-#include "EventDebugger.hh"
-
-// slic
-#include "EventDebuggerMessenger.hh"
-
-// geant4
-#include "G4Run.hh"
-#include "G4Event.hh"
-#include "G4UImanager.hh"
-
-namespace slic {
-EventDebugger::EventDebugger() :
- Module("EventDebugger"), m_debugging(false), m_haveDebugMacro(false), m_haveCleanupMacro(false), m_forceDebugMode(
- false) {
- m_messenger = new EventDebuggerMessenger();
-}
-
-EventDebugger::~EventDebugger() {
- delete m_messenger;
-}
-
-void EventDebugger::enableDebug(bool e) {
- m_forceDebugMode = e;
-}
-
-void EventDebugger::addDebugEvent(G4int eventNum) {
- if (!haveDebugEvent(eventNum)) {
- m_events.push_back(eventNum);
-#ifdef SLIC_LOG
- log() << LOG::okay << "EventDebugger - Added debug event <" << eventNum << ">" << LOG::done;
-#endif
- }
-#ifdef SLIC_LOG
- else {
- log() << LOG::warning << "EventDebugger - Ignoring dup event #" << eventNum << LOG::done;
- }
-#endif
-}
-
-void EventDebugger::setDebugMacro(std::string mac) {
- m_haveDebugMacro = true;
- m_debugMacro = mac;
-}
-
-void EventDebugger::setCleanupMacro(std::string mac) {
- m_haveCleanupMacro = true;
- m_cleanupMacro = mac;
-}
-
-void EventDebugger::clearDebugEvents() {
- m_events.clear();
-}
-
-void EventDebugger::sortDebugEvents() {
- std::sort(m_events.begin(), m_events.end());
-}
-
-bool EventDebugger::haveDebugEvent(G4int evtNum) const {
- for (DebugEventList::const_iterator it = m_events.begin(); it != m_events.end(); it++) {
- if ((*it) == evtNum) {
- return true;
- }
- }
- return false;
-}
-
-void EventDebugger::beginRun(const G4Run*) {
- sortDebugEvents();
-}
-
-void EventDebugger::beginEvent(const G4Event* evt) {
- if (m_forceDebugMode || haveDebugEvent(evt->GetEventID())) {
- m_debugging = true;
- execDebugMacro();
- }
-}
-
-void EventDebugger::endEvent(const G4Event*) {
- if (m_debugging) {
- execCleanupMacro();
- m_debugging = false;
- }
-}
-
-void EventDebugger::execDebugMacro() {
- if (m_haveDebugMacro) {
- G4UImanager::GetUIpointer()->ApplyCommand("/control/execute " + m_debugMacro);
-
-#ifdef SLIC_LOG
- if ( !m_haveCleanupMacro ) {
- log() << "WARNING: Executed debugging macro, but cleanup macro was not set." << LOG::done;
- }
-#endif
- }
-#ifdef SLIC_LOG
- else {
- log() << LOG::error << "ERROR: debug macro not set" << LOG::done;
- }
-#endif
-}
-
-void EventDebugger::execCleanupMacro() {
- if (m_haveCleanupMacro) {
- G4UImanager::GetUIpointer()->ApplyCommand("/control/execute " + m_cleanupMacro);
- }
-#ifdef SLIC_LOG
- else {
- log() << LOG::error << "ERROR: cleanup macro not set" << LOG::done;
- }
-#endif
-}
-} // namespace slic
projects/slic/trunk/src
--- projects/slic/trunk/src/EventDebuggerMessenger.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/EventDebuggerMessenger.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,68 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/EventDebuggerMessenger.cc,v 1.6 2013-06-26 01:57:02 jeremy Exp $
-#include "EventDebuggerMessenger.hh"
-
-// slic
-#include "EventDebugger.hh"
-
-// lcdd
-#include "lcdd/util/StringUtil.hh"
-
-// geant4
-#include "G4UIcmdWithAString.hh"
-#include "G4UIcmdWithABool.hh"
-
-namespace slic {
-EventDebuggerMessenger::EventDebuggerMessenger() {
- defineCommands();
-}
-
-EventDebuggerMessenger::~EventDebuggerMessenger() {
-}
-
-void EventDebuggerMessenger::SetNewValue(G4UIcommand* cmd, G4String newVals) {
- EventDebugger* debugger = EventDebugger::instance();
- if (cmd == m_debugMacroCmd) {
- debugger->setDebugMacro(newVals);
- } else if (cmd == m_cleanupMacroCmd) {
- debugger->setCleanupMacro(newVals);
- } else if (cmd == m_addDebugEventsCmd) {
-
- std::vector<std::string> strList;
- const std::string str = std::string(newVals);
- const std::string delim = " ";
-
- StringUtil::split(str, delim, strList);
-
- EventDebugger::DebugEventList dbgList;
- for (std::vector<std::string>::iterator it = strList.begin(); it != strList.end(); it++) {
- debugger->addDebugEvent(StringUtil::toInt(*it));
- }
- } else if (cmd == m_clearDebugEventsCmd) {
- debugger->clearDebugEvents();
- } else if (cmd == m_enableDebugCmd) {
- debugger->enableDebug(m_enableDebugCmd->GetNewBoolValue(newVals.c_str()));
- }
-}
-
-void EventDebuggerMessenger::defineCommands() {
- m_debugDir = new G4UIdirectory("/debug/");
- m_debugDir->SetGuidance("Debugging commands. [SLIC]");
-
- m_debugMacroCmd = new G4UIcmdWithAString("/debug/debugMacro", this);
- m_debugMacroCmd->SetGuidance("Set macro for debugging.");
-
- m_cleanupMacroCmd = new G4UIcmdWithAString("/debug/cleanupMacro", this);
- m_cleanupMacroCmd->SetGuidance("Set macro to cleanup debug state.");
-
- m_addDebugEventsCmd = new G4UIcmdWithAString("/debug/addEvents", this);
- m_addDebugEventsCmd->SetGuidance("Add event IDs to debug, separated by spaces.");
-
- m_clearDebugEventsCmd = new G4UIcommand("/debug/clearEvents", this);
- m_clearDebugEventsCmd->SetGuidance("Clear the list of events to debug.");
-
- m_enableDebugCmd = new G4UIcmdWithABool("/debug/enableDebug", this);
- m_enableDebugCmd->SetGuidance("Enable event debugging for all subsequent events.");
- m_enableDebugCmd->SetParameterName("enable", true);
- m_enableDebugCmd->SetDefaultValue(true);
-}
-} // namespace slic
projects/slic/trunk/src
--- projects/slic/trunk/src/EventSource.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/EventSource.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -11,7 +11,7 @@
}
EventSource::EventSource(const std::string& name) :
- Module(name), m_name(name) {
+ Module(name), m_name(name), nEventsGenerated(0) {
}
void EventSource::reset() {
projects/slic/trunk/src
--- projects/slic/trunk/src/EventSourceManager.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/EventSourceManager.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,31 +1,33 @@
// $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 "EventSourceManager.hh"
-#include "GPSEventSource.hh"
-#include "ParticleGunEventSource.hh"
-#include "StdHepEventSource.hh"
-#include "LcioEventSource.hh"
#include "FileUtil.hh"
+
#include "EventAction.hh"
+#include "RunManager.hh"
#include "SlicApplication.hh"
-#include "RunManager.hh"
+#include "GPSEventSource.hh"
+#include "LcioEventSource.hh"
+#include "ParticleGunEventSource.hh"
+#include "StdHepEventSource.hh"
+
// LCDD
#include "lcdd/util/StringUtil.hh"
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), m_newFilename(
- false), m_ngen(0), m_nskip(0), m_newSource(true), m_sourceType(eUnknown), m_lorentzTransformationAngle(
- 0.), m_zSmearingParam(0.) {
+ Module("EventSourceManager"), m_currentEventSource(0), m_filename(""), m_fileIsSet(false),
+ m_newFilename(false), m_ngen(0), m_nskip(0), m_newSource(true), m_sourceType(eUnknown),
+ m_lorentzTransformationAngle(0.), m_zSmearingParam(0.), _enablePrintEvent(false) {
// messenger with generator command macro bindings
m_messenger = new GeneratorMessenger();
@@ -98,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 {
@@ -115,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;
}
@@ -137,7 +139,7 @@
}
void EventSourceManager::dumpCurrentEvent() {
- m_currentEventSource->dumpCurrentEvent();
+ m_currentEventSource->printCurrentEvent();
}
void EventSourceManager::printNumEventsGenerated() {
@@ -241,6 +243,10 @@
void EventSourceManager::endEvent(const G4Event* anEvent) {
m_currentEventSource->endEvent(anEvent);
+
+ if (_enablePrintEvent) {
+ m_currentEventSource->printCurrentEvent();
+ }
}
void EventSourceManager::beginRun(const G4Run* aRun) {
@@ -261,10 +267,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 +292,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 +341,7 @@
return gun;
}
-void EventSourceManager::GeneratePrimaryVertex(G4Event* evt) {
+void EventSourceManager::generateNextEvent(G4Event* evt) {
generate(evt);
}
@@ -346,10 +349,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/trunk/src
--- projects/slic/trunk/src/G4MagIntegratorStepperFactory.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/G4MagIntegratorStepperFactory.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -2,6 +2,7 @@
#include "G4MagIntegratorStepperFactory.hh"
// G4
+#include "G4SystemOfUnits.hh"
#include "G4Mag_UsualEqRhs.hh"
#include "G4ChordFinder.hh"
#include "G4MagIntegratorStepper.hh"
projects/slic/trunk/src
--- projects/slic/trunk/src/GPSEventSource.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/GPSEventSource.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -35,7 +35,7 @@
}
/** This is mostly wrong until after the first event generation. */
-void GPSEventSource::dumpCurrentEvent() {
+void GPSEventSource::printCurrentEvent() {
log() << LOG::okay << "Dumping GPS Event info..." << LOG::done;
G4ParticleDefinition* pdef = m_generator->GetParticleDefinition();
projects/slic/trunk/src
--- projects/slic/trunk/src/GeneratorMessenger.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/GeneratorMessenger.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -58,8 +58,9 @@
mgr->setFilename(s);
}
// dumpCurrentEvent
- else if (cmd == m_dumpCurrentEventCmd) {
- mgr->dumpCurrentEvent();
+ else if (cmd == m_enablePrintEventCmd) {
+ G4bool enablePrintEvent = StringUtil::toBool(s);
+ mgr->enablePrintEvent(true);
}
// printNumEventsGenerated
else if (cmd == m_printNumEventsGeneratedCmd) {
@@ -124,8 +125,8 @@
m_filenameCmd->SetParameter(p);
// dumpCurrentEvent
- m_dumpCurrentEventCmd = new G4UIcommand("/generator/dumpEvent", this);
- m_dumpCurrentEventCmd->SetGuidance("Dump information about the current generator event.");
+ m_enablePrintEventCmd = new G4UIcommand("/generator/printEvent", this);
+ m_enablePrintEventCmd->SetGuidance("Dump information about the current generator event.");
// printNumEventsGenerated
m_printNumEventsGeneratedCmd = new G4UIcommand("/generator/printNumEventsGenerated", this);
projects/slic/trunk/src
--- projects/slic/trunk/src/LCExtendedParticles.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LCExtendedParticles.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -6,7 +6,6 @@
// geant4
#include "G4ProcessManager.hh"
#include "G4SystemOfUnits.hh"
-#include "G4ParticleTableIterator.hh"
namespace slic {
@@ -69,12 +68,9 @@
}
void LCExtendedParticles::ConstructProcess() {
-
- G4ParticleTable::G4PTblDicIterator* particleIterator = G4ParticleTable::GetParticleTable()->GetIterator();
- particleIterator->reset();
-
- while ((*particleIterator)()) {
- G4ParticleDefinition* pdef = particleIterator->value();
+ aParticleIterator->reset();
+ while ((*aParticleIterator)()) {
+ G4ParticleDefinition* pdef = aParticleIterator->value();
G4ProcessManager* pmgr = pdef->GetProcessManager();
if (pdef->GetParticleType() == "extended" || pdef->GetParticleType() == "susy") {
projects/slic/trunk/src
--- projects/slic/trunk/src/LCOpticalPhysics.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LCOpticalPhysics.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -55,12 +55,11 @@
//G4OpticalSurfaceModel themodel = unified;
//theBoundaryProcess->SetModel(themodel);
- G4ParticleTable::G4PTblDicIterator* particleIterator = G4ParticleTable::GetParticleTable()->GetIterator();
- particleIterator->reset();
+ aParticleIterator->reset();
- while ((*particleIterator)()) {
+ while ((*aParticleIterator)()) {
- G4ParticleDefinition* particle = particleIterator->value();
+ G4ParticleDefinition* particle = aParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
projects/slic/trunk/src
--- projects/slic/trunk/src/LCSUSYPhysics.cc (rev 0)
+++ projects/slic/trunk/src/LCSUSYPhysics.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -0,0 +1,75 @@
+#include "LCSUSYPhysics.hh"
+
+// slic
+#include "HepPDTManager.hh"
+#include "LCSUSYDecay.hh"
+
+// geant4
+#include "G4ProcessManager.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4UnknownDecay.hh"
+#include "G4SystemOfUnits.hh"
+
+namespace slic {
+LCSUSYPhysics::LCSUSYPhysics(const G4String& name) :
+ G4VPhysicsConstructor(name) {
+ ;
+}
+
+LCSUSYPhysics::~LCSUSYPhysics() {
+ ;
+}
+
+void LCSUSYPhysics::ConstructParticle() {
+ // Get the particle data table to be imported.
+ HepPDT::ParticleDataTable* pdtbl = HepPDTManager::instance()->getParticleDataTable();
+
+ // Loop over all particles.
+ for (HepPDT::ParticleDataTable::const_iterator it = pdtbl->begin(); it != pdtbl->end(); it++) {
+ HepPDT::ParticleID id = it->first;
+ HepPDT::ParticleData pdata = it->second;
+
+ // Only adding SUSY particles.
+ if (id.isSUSY()) {
+ // Create a G4ParticleDefinition for this SUSY particle.
+ new G4ParticleDefinition(pdata.PDTname(), // name
+ pdata.mass().value() * GeV, // mass
+ pdata.totalWidth().value(), // width
+ pdata.charge(), // charge
+ 0, // 2*spin
+ 0, // parity
+ 0, // C-conjugation
+ 0, // 2*isospin
+ 0, // 2*isospin3
+ 0, // G-parity
+ "susy", // type
+ 0, // lepton number
+ 0, // baryon number
+ pdata.pid(), // PDG encoding
+ false, // stable
+ 1e30, // lifetime (DEBUG)
+ 0, // decay table
+ false // short lived
+ );
+ }
+ }
+}
+
+void LCSUSYPhysics::ConstructProcess() {
+ aParticleIterator->reset();
+ while ((*aParticleIterator)()) {
+ G4ParticleDefinition* pdef = aParticleIterator->value();
+ G4ProcessManager* pmgr = pdef->GetProcessManager();
+ if (m_decay.IsApplicable(*pdef)) {
+ if (pdef->GetPDGCharge() != 0) {
+ pmgr->AddProcess(&m_scatter, -1, 1, 1); // multiple scattering
+ pmgr->AddProcess(&m_ionise, -1, 2, 2); // ionisation
+ pmgr->AddProcess(&m_decay, -1, -1, 3); // decay
+ } else {
+ pmgr->AddProcess(&m_scatter, -1, 1, 1); // multiple scattering
+ pmgr->AddProcess(&m_decay, -1, -1, 2); // decay
+ }
+ }
+ }
+}
+}
projects/slic/trunk/src
--- projects/slic/trunk/src/LcioEventSource.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LcioEventSource.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,93 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/LcioEventSource.cc,v 1.10 2013-11-06 00:23:37 jeremy Exp $
-#include "LcioEventSource.hh"
-
-// slic
-#include "LcioMcpFilter.hh"
-
-namespace slic {
-
-LcioEventSource::LcioEventSource(const std::string& fname) :
- EventSourceWithInputFile("LcioEventSource", fname), m_event(0) {
- m_reader = IOIMPL::LCFactory::getInstance()->createLCReader();
-
- m_filter = new LcioMcpFilter();
-}
-
-LcioEventSource::~LcioEventSource() {
- delete m_reader;
-}
-
-// open the current file
-void LcioEventSource::open() {
- try {
-#ifdef SLIC_LOG
- log() << LOG::verbose << "opening filename <" << m_filename << ">" << LOG::done;
-#endif
- m_reader->open(m_filename);
- m_fileIsOpen = true;
- } catch (IO::IOException& ioe) {
- log() << LOG::error << ioe.what() << LOG::done;
- G4Exception("", "", FatalException, "Caught exception while reading LCIO input file.");
- }
-}
-
-// close the current file
-void LcioEventSource::close() {
- try {
- m_reader->close();
- m_fileIsOpen = false;
- } catch (...) {
- }
-}
-
-// read the next event
-void LcioEventSource::readNextEvent() {
- m_event = m_reader->readNextEvent();
-
- if (0 == m_event) {
- m_eof = true;
- }
-}
-
-void LcioEventSource::dumpCurrentEvent() {
- log() << LOG::debug << "LcioEventSource - Initial MCParticle Collection" << LOG::done;
-
- LcioMcpManager::instance()->printMcpCollection("initial");
-}
-
-void LcioEventSource::generate(G4Event* anEvent) {
- assert( LcioMcpManager::instance()->getInitialMcpCollection());
- LcioManager::instance()->getGenerator()->generatePrimaryVertexFromMcpCollection(
- LcioMcpManager::instance()->getInitialMcpCollection(), 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()->findParticleCollection(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/trunk/src
--- projects/slic/trunk/src/LcioFileNamer.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LcioFileNamer.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -9,6 +9,7 @@
#include "RunManager.hh"
#include "TimeUtil.hh"
#include "PhysicsListManager.hh"
+#include "Geant4VersionInfo.hh"
// LCDD
#include "lcdd/core/LCDDProcessor.hh"
@@ -43,8 +44,10 @@
std::string LcioFileNamer::getFieldValue(std::string field) {
std::string value;
+ static std::string sep("-");
+
if (field == "application" || field == "app") {
- value = PackageInfo::getAbbrevString() + "-" + PackageInfo::getVersionString();
+ value = PackageInfo::getShortName() + "-v" + PackageInfo::getVersion(sep);
} else if (field == "geometry" || field == "geo") {
value = LCDDProcessor::instance()->getDetectorName();
} else if (field == "date") {
@@ -60,7 +63,7 @@
} else if (field == "physics") {
value = PhysicsListManager::instance()->getCurrentListName();
} else if (field == "geant4") {
- value = SlicApplication::instance()->getGeant4VersionString();
+ value = "Geant4-v" + Geant4VersionInfo::getVersion(sep);
} else {
#ifdef SLIC_LOG
log() << LOG::error << "Ignoring unknown autoname field <" << field << ">." << LOG::done;
projects/slic/trunk/src
--- projects/slic/trunk/src/LcioHitsCollectionBuilder.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LcioHitsCollectionBuilder.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -2,8 +2,8 @@
// SLIC
#include "LcioHitsCollectionBuilder.hh"
-#include "LcioMcpManager.hh"
#include "HitsCollectionUtil.hh"
+#include "MCParticleManager.hh"
#include "SlicApplication.hh"
// LCDD
@@ -14,18 +14,18 @@
#include "G4SystemOfUnits.hh"
// LCIO
+#include "EVENT/MCParticle.h"
#include "IMPL/LCCollectionVec.h"
#include "IMPL/LCFlagImpl.h"
-#include "IMPL/MCParticleImpl.h"
// LCIO
using IMPL::SimCalorimeterHitImpl;
using IMPL::SimTrackerHitImpl;
-using IMPL::MCParticleImpl;
using IMPL::LCEventImpl;
using IMPL::LCFlagImpl;
using IMPL::LCCollectionVec;
using EVENT::LCIO;
+using EVENT::MCParticle;
// STL
#include <iostream>
@@ -36,8 +36,6 @@
LcioHitsCollectionBuilder::LcioHitsCollectionBuilder() :
Module("LcioHitsCollectionBuilder"), m_storeMomentum(0) {
- // set local LcioMcpManager ptr
- m_mcpManager = LcioMcpManager::instance();
// setup default coll flag for cal hits
setFlagDefaults();
@@ -288,20 +286,14 @@
simTrackerHit->setTime(tEdep);
// Cell ID.
-//#if LCIO_VERSION_GE(1, 60)
- // New method for 64-bit IDs.
simTrackerHit->setCellID0(trackerHit->getId());
-//#else
- // Old method for 32-bit IDs.
- //simTrkHit->setCellID(trkHit->getId());
-//#endif
// MCP using McpManager
- MCParticleImpl* mcp = m_mcpManager->getMaps()->findParticle(trackerHit->getTrackID());
+ MCParticle* mcp = MCParticleManager::instance()->findMCParticle(trackerHit->getTrackID());
if (!mcp) {
log().error(
- "No MCP found for trackID <" + StringUtil::toString(trackerHit->getTrackID()) + "> for trk hit.");
+ "No MCParticle found for trackID <" + StringUtil::toString(trackerHit->getTrackID()) + "> when processing SimTrackerHit");
} else {
simTrackerHit->setMCParticle(mcp);
}
@@ -332,37 +324,29 @@
const HitContribution& contrib = (*iter);
// Get the MCParticle pointer from the track ID.
- MCParticleImpl* contribMcp = m_mcpManager->getMaps()->findParticle(contrib.getTrackID());
+ MCParticle* contribMcp = MCParticleManager::instance()->findMCParticle(contrib.getTrackID());
if (contribMcp != 0) {
// Add the MCParticle contribution to the hit.
-//#if LCIO_VERSION_GE(1, 60)
- // Newer LCIO versions have the step position for contributions.
simCalHit->addMCParticleContribution(
contribMcp,
(float)(contrib.getEdep()/GeV),
(float)(contrib.getGlobalTime()),
contrib.getPDGID(),
const_cast<float*>(contrib.getPosition()));
-//#else
- // Older LCIO versions do not include the step position.
-// simCalHit->addMCParticleContribution(contribMcp, (float) (contrib.getEdep() / GeV),
-// (float) contrib.getGlobalTime(), contrib.getPDGID());
-//#endif
++ncontrib;
}
// Problem! Contributing particle is missing from MCParticle list.
#ifdef SLIC_LOG
else
{
- log() << LOG::always << "ERROR: Could not find MCParticle from track ID <" << contrib.getTrackID() << ">." << LOG::endl;
+ log() << LOG::error << "No MCParticle from track ID <" << contrib.getTrackID() << "> when processing SimCalorimeterHit" << LOG::endl;
}
#endif
}
#ifdef SLIC_LOG
- if ( ncontrib == 0 )
- {
+ if ( ncontrib == 0 ) {
log().error("No hit contributions for CalorimeterHit.");
}
#endif
projects/slic/trunk/src
--- projects/slic/trunk/src/LcioManager.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LcioManager.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,26 +1,25 @@
// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/LcioManager.cc,v 1.96 2013-11-06 00:23:37 jeremy Exp $
#include "LcioManager.hh"
-// slic
+// SLIC
#include "EventSourceManager.hh"
#include "StdHepEventSource.hh"
#include "FileUtil.hh"
#include "TimeUtil.hh"
#include "PackageInfo.hh"
#include "LcioMessenger.hh"
-#include "LcioMcpManager.hh"
#include "LcioFileNamer.hh"
-#include "TrajectoryManager.hh"
#include "SlicApplication.hh"
#include "RunManager.hh"
+#include "Geant4VersionInfo.hh"
-// lcdd
+// LCDD
#include "lcdd/hits/CalorimeterHit.hh"
#include "lcdd/hits/TrackerHit.hh"
#include "lcdd/core/LCDDProcessor.hh"
#include "lcdd/util/StringUtil.hh"
-// lcio
+// LCIO
#include "EVENT/LCIO.h"
#include "IMPL/LCCollectionVec.h"
#include "IOIMPL/LCFactory.h"
@@ -28,7 +27,7 @@
#include "UTIL/LCTOOLS.h"
#include "UTIL/lStdHep.hh"
-// geant4
+// Geant4
#include "G4EventManager.hh"
#include "G4Run.hh"
#include "G4RunManager.hh"
@@ -36,7 +35,7 @@
#include "G4VHitsCollection.hh"
#include "G4SDManager.hh"
-// std
+// STL
#include <ctime>
using IMPL::LCRunHeaderImpl;
@@ -51,21 +50,14 @@
namespace slic {
-string LcioManager::m_defaultFileName = "outfile";
+string LcioManager::m_defaultFileName = "slicEvents";
LcioManager::LcioManager() :
- Module("LcioManager"), m_McpFinalColl(0), m_writer(0), m_runHdr(0), m_fileExistsAction(
- LcioManager::eFail), m_runNumber(0), m_enableDumpEvent(false), m_writerIsOpen(false), m_usingAutoname(
- false) {
+ Module("LcioManager"), m_writer(0), m_runHdr(0), m_fileExistsAction(LcioManager::eFail),
+ m_runNumber(0), m_enableDumpEvent(false), m_writerIsOpen(false), m_usingAutoname(false) {
// Initialize the Geant4 UI messenger for the LCIO.
m_messenger = new LcioMessenger(this);
- // Initialize the LCIO-based event generator.
- m_eventGenerator = new LcioPrimaryGenerator(this);
-
- // Initialize the LCIO MCParticle manager.
- m_mcpManager = LcioMcpManager::instance();
-
// Initialize the LCIO HitsCollectionBuilder.
m_HCBuilder = new LcioHitsCollectionBuilder();
@@ -80,12 +72,12 @@
}
LcioManager::~LcioManager() {
- deleteWriter();
+ deleteWriter();
+
if (m_messenger != 0) {
delete m_messenger;
}
-
}
void LcioManager::openLcioFile() {
@@ -251,11 +243,10 @@
m_runHdr = new LCRunHeaderImpl();
// Write SLIC version into run header.
- m_runHdr->parameters().setValue("SLIC_VERSION", PackageInfo::getVersionString());
+ m_runHdr->parameters().setValue("SLIC_VERSION", PackageInfo::getVersion());
// Write Geant4 version into run header.
- m_runHdr->parameters().setValue("GEANT4_VERSION",
- SlicApplication::instance()->getGeant4VersionString().replace(0, 7, ""));
+ m_runHdr->parameters().setValue("GEANT4_VERSION", Geant4VersionInfo::getVersion());
// set run number
m_runHdr->setRunNumber(m_runNumber);
@@ -329,43 +320,42 @@
}
}
-const string& LcioManager::getPath() const {
- return m_path;
-}
-
-const string& LcioManager::getFilename() const {
- return m_filename;
-}
-
LCEventImpl* LcioManager::createLCEvent(const G4Event* anEvent) {
- assert( anEvent);
+ /* Create the LCEvent and set basic parameters. */
LCEventImpl* lcevt = new LCEventImpl();
lcevt->setEventNumber(anEvent->GetEventID());
lcevt->setRunNumber(m_runHdr->getRunNumber());
lcevt->setDetectorName(m_runHdr->getDetectorName());
- // 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);
- }
+ /* Get the current MCParticle collection. */
+ LCCollection* mcpColl = MCParticleManager::instance()->getMCParticleCollection();
- // Write Geant4 version into event header.
- lcevt->parameters().setValue("GEANT4_VERSION",
- SlicApplication::instance()->getGeant4VersionString().replace(0, 7, ""));
+ if (mcpColl != 0) {
- // Write SLIC version into event header.
- lcevt->parameters().setValue("SLIC_VERSION", PackageInfo::getVersionString());
+ /* Set the event weight. */
+ float eventWeight = mcpColl->getParameters().getFloatVal("_weight");
+ if (eventWeight != 0.0) {
+ lcevt->setWeight(eventWeight);
+ }
+ /* Set the idrup. */
+ int idrup = mcpColl->getParameters().getIntVal("_idrup");
+ //G4cout << "idrup = " << idrup << G4endl;
+ if (idrup != 0) {
+ lcevt->parameters().setValue("_idrup", idrup);
+ }
+ }
+
+ // Write SLIC version into event header.
+ lcevt->parameters().setValue("SLIC_VERSION", PackageInfo::getVersion());
+
+ // Write Geant4 version into event header.
+ lcevt->parameters().setValue("GEANT4_VERSION", Geant4VersionInfo::getVersion());
+
+ /* Set the current LCEvent. */
setCurrentLCEvent(lcevt);
+
return lcevt;
}
@@ -374,43 +364,27 @@
}
void LcioManager::endEvent(const G4Event*) {
- if (!RunManager::instance()->isRunAborted()) {
- // create LCEvent
- createLCEvent();
+ // create HC in current LCEvent from current G4Event using builder
+ createHitsCollections();
- // create Mcp coll in LcioMcpManager
- m_mcpManager->endEvent(G4EventManager::GetEventManager()->GetNonconstCurrentEvent());
+ // Dump event stats.
+ if (m_enableDumpEvent) {
+ LCTOOLS::dumpEventDetailed(m_currentLCEvent);
+ }
- // create mcp collection from input event
- createFinalMcpCollection();
+ // set timestamp
+ setEventTimeStamp();
- // If selected, add the initial MCParticle collection to the event.
- if (LcioMcpManager::instance()->writeInitialMCParticleCollection()) {
- addInitialMCParticleCollection();
- }
+ // write event
+ m_writer->writeEvent(m_currentLCEvent);
- // create HC in current LCEvent from current G4Event using builder
- createHitsCollections();
+ // flush writer
+ m_writer->flush();
- // Dump event stats.
- if (m_enableDumpEvent) {
- LCTOOLS::dumpEventDetailed(m_currentLCEvent);
- //LCTOOLS::printMCParticles( LcioMcpManager::instance()->getFinalMcpCollection() );
- }
-
- // set timestamp
- setEventTimeStamp();
-
- // write event
- m_writer->writeEvent(m_currentLCEvent);
-
- // flush writer
- m_writer->flush();
-
- // delete event's transient objects, including current LCEvent
- reset();
- }
+ // delete event's transient objects, including current LCEvent
+ delete m_currentLCEvent;
+ m_currentLCEvent = 0;
}
void LcioManager::setEventTimeStamp() {
@@ -421,20 +395,6 @@
m_HCBuilder->createHitCollectionsFromEvent(G4EventManager::GetEventManager()->GetNonconstCurrentEvent(), m_currentLCEvent);
}
-void LcioManager::reset() {
- // delete transient event container
- delete m_currentLCEvent;
- m_currentLCEvent = 0;
-
- // reset mcp mgr, including clearing maps
- m_mcpManager->reset();
-}
-
-void LcioManager::createFinalMcpCollection() {
- // add Mcp coll to current event
- getCurrentLCEvent()->addCollection(m_mcpManager->getFinalMcpCollection(), LCIO::MCPARTICLE);
-}
-
void LcioManager::addCollection(EVENT::LCEvent* event, EVENT::LCCollection* collection,
const std::string& collectionName) {
event->addCollection(collection, collectionName);
@@ -444,8 +404,4 @@
getCurrentLCEvent()->addCollection(collection, collectionName);
}
-void LcioManager::addInitialMCParticleCollection() {
- const std::string& name = std::string(LCIO::MCPARTICLE) + std::string("Initial");
- addCollection(LcioMcpManager::instance()->getInitialMcpCollection(), name);
}
-}
projects/slic/trunk/src
--- projects/slic/trunk/src/LcioMcpFactory.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LcioMcpFactory.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,659 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/LcioMcpFactory.cc,v 1.25 2013-11-06 00:23:37 jeremy Exp $
-#include "LcioMcpFactory.hh"
-
-// slic
-#include "EventSourceManager.hh"
-#include "LcioManager.hh"
-#include "LcioMcpManager.hh"
-#include "LcioMcpStatusSetter.hh"
-#include "LcioMcpUtil.hh"
-#include "Trajectory.hh"
-#include "TrajectoryManager.hh"
-
-// geant4
-#include "globals.hh"
-#include "G4Event.hh"
-#include "G4EventManager.hh"
-#include "G4TrajectoryContainer.hh"
-#include "G4SystemOfUnits.hh"
-
-// lcio
-#include "IMPL/LCCollectionVec.h"
-#include "IMPL/LCGenericObjectImpl.h"
-
-using IMPL::LCCollectionVec;
-using IMPL::LCGenericObjectImpl;
-using IMPL::MCParticleImpl;
-using EVENT::LCIO;
-using EVENT::MCParticle;
-
-namespace slic {
-LcioMcpFactory::LcioMcpFactory(LcioMcpManager* manager) :
- Module("LcioMcpFactory"), m_manager(manager), m_finalColl(0), m_currentTrajectoryContainer(
- 0) {
-}
-
-LcioMcpFactory::~LcioMcpFactory() {
-}
-
-void LcioMcpFactory::createFinalMcpCollection(const G4Event* event) {
- // new coll for final Mcps
- m_finalColl =
- static_cast<LCCollectionVec*>(m_manager->getFinalMcpCollection());
-
- // Set current G4TrajectoryContainer
- m_currentTrajectoryContainer = event->GetTrajectoryContainer();
-
- // Check if using StdHep or LCIO e.g. a file-based source?
- if (EventSourceManager::instance()->isFileSource()) {
- // Create from the initial coll of Mcp.
- createFinalCollectionFromInitial(
- m_manager->getInitialMcpCollection());
- }
- // G4 source using trajectories only.
- else {
- // Create from trajectory container of current G4Event.
- createFinalCollectionFromTrajectories(
- m_currentTrajectoryContainer);
- }
-
- // Fill supplementary collection with MCParticle endpoint energies.
- fillEndPointEnergy(m_finalColl);
-}
-
-IMPL::MCParticleImpl* LcioMcpFactory::fromTrajectory(Trajectory* trj) {
-#ifdef SLIC_LOG
- log() << LOG::debug << "createMcpFromTrajectory() - trackId: " << trj->GetTrackID() << LOG::done;
- log() << LOG::debug << "trj PDG: " << trj->GetParticleDefinition()->GetPDGEncoding() << LOG::done;
-#endif
-
- // new Mcp
- MCParticleImpl* mcp = new MCParticleImpl();
-
- // PDG
- mcp->setPDG(trj->GetPDGEncoding());
-
- // mom
- G4ThreeVector mom = trj->GetInitialMomentum() / GeV;
- float fmom[3] = { mom.x(), mom.y(), mom.z() };
- mcp->setMomentum(fmom);
-
- // Get particle def.
- const G4ParticleDefinition* pdef = trj->GetParticleDefinition();
-
- // To store mass for new particle.
- G4double m = 0;
-
- /*
- * For particles from extended particle table supplied by slic, use the StdHep mass instead.
- * Also, use the StdHep mass for PDG codes Geant4 doesn't know that aren't in particle.tbl.
- */
- if (pdef == 0
- || (pdef->GetParticleType().compare("extended") == 0
- || pdef->GetParticleType().compare("susy") == 0)) {
- // Get the primary particle from the LCIO data map.
- MCParticle* primaryMcp =
- LcioMcpManager::instance()->getMaps()->findInitialParticle(
- trj);
-
- // Require that this particle have a matching primary MCParticle to get mass.
- if (primaryMcp != 0) {
- m = primaryMcp->getMass();
-#ifdef SLIC_LOG
- log() << "mass from primary: " << m << LOG::endl;
-#endif
- }
- // Need a primary particle to get mass.
- else {
- G4Exception("", "", FatalException,
- "Could not find required LCIO primary particle.");
- }
- }
- // Use Geant4's mass from the particle definition, which is appropriate for SM particles.
- else {
- m = trj->GetParticleDefinition()->GetPDGMass() / GeV;
-#ifdef SLIC_LOG
- log() << "mass from pdef: " << m << LOG::endl;
-#endif
- }
-
- // mass
- mcp->setMass(m);
-
- // charge
- mcp->setCharge(trj->GetCharge());
-
- // begin and end points
- G4int npts = trj->GetPointEntries();
-
- // begin pnt
- G4VTrajectoryPoint* beginTrjPnt = trj->GetPoint(0);
- G4ThreeVector beginPos = beginTrjPnt->GetPosition() / mm;
- double fBeginPos[3] = { beginPos.x(), beginPos.y(), beginPos.z() };
- mcp->setVertex(fBeginPos);
-
- // end pnt
- G4VTrajectoryPoint* endTrjPnt = trj->GetPoint(npts - 1);
- G4ThreeVector endPos = endTrjPnt->GetPosition() / mm;
- double fEndPos[3] = { endPos.x(), endPos.y(), endPos.z() };
- mcp->setEndpoint(fEndPos);
-
- // time
- mcp->setTime(trj->GetGlobalTime());
-
- // set status codes from trajectory info (with no dep on any associated Mcps)
- LcioMcpStatusSetter::setStatusCodesFromTrajectory(trj, mcp);
-
- // add to map of trkID -> mcp
-#ifdef SLIC_LOG
- log() << LOG::debug << "adding trackId <" << trj->GetTrackID() << "> to MCP <" << mcp << "> link." << LOG::done;
-#endif
- m_manager->getMaps()->addTrackToMCParticle(trj->GetTrackID(), mcp);
-
- return mcp;
-}
-
-// FIXME: This method doesn't seem to work correctly!
-void LcioMcpFactory::createFinalCollectionFromTrajectories(
- G4TrajectoryContainer* m_currentTrajectoryContainer) {
- if (m_currentTrajectoryContainer) {
- // loop over trj cont
- int n_trj = m_currentTrajectoryContainer->entries();
- for (int i = 0; i < n_trj; i++) {
- Trajectory* trj =
- static_cast<Trajectory*>((*m_currentTrajectoryContainer)[i]);
-
- // Only add primary trajectories; daughters are picked up within.
- if (trj->GetParentID() == 0) {
-#ifdef SLIC_LOG
- log() << LOG::debug << "Making primary MCParticle for trajectory with trackID <" << trj->GetTrackID() << ">." << LOG::done;
-#endif
-
- // Create MCParticle output from trajectory.
- MCParticleImpl* mcp = fromTrajectory(trj);
-
- // Add to the final MCParticle collection.
- m_finalColl->addElement(mcp);
-
- // Process the daughters.
-#ifdef SLIC_LOG
- log() << LOG::debug << "Making primary MCParticle for trajectory with trackID <" << trj->GetTrackID() << ">." << LOG::done;
-#endif
- createDaughters(mcp, trj->GetTrackID());
- }
- }
- } else {
- G4Exception("", "", FatalException, "G4TrajectoryContainer is null!");
- }
-
- // DEBUG
- LcioMcpManager::instance()->getMaps()->printTracksToParticles();
-
-}
-
-void LcioMcpFactory::createDaughters(
- MCParticleImpl* parMcp, int parTrkID) {
-#ifdef SLIC_LOG
- log() << LOG::debug << "addMcpDaughtersFromTraj - parTrkId <" << parTrkID << ">." << LOG::done;
-#endif
-
- // loop over trj cont
- int n_trj = m_currentTrajectoryContainer->entries();
- for (int i = 0; i < n_trj; i++) {
- Trajectory* trj =
- static_cast<Trajectory*>((*m_currentTrajectoryContainer)[i]);
-
- if (trj->GetParentID() == parTrkID) {
-
- // lkp dau Mcp by track ID
- MCParticleImpl* dauMcp = m_manager->getMaps()->findParticle(
- trj->GetTrackID());
-
- if (dauMcp == 0) {
- dauMcp = fromTrajectory(trj);
- }
-
- if (dauMcp == 0)
- G4Exception("", "", FatalException,
- "Failed to create MCParticle.");
-
- // add this to final coll
- m_finalColl->addElement(dauMcp);
-
- // Recursively add daughters of this daughter.
- createDaughters(dauMcp, trj->GetTrackID());
-
- // Set the parent particle in the daughter.
- dauMcp->addParent(parMcp);
- }
- }
-}
-
-void LcioMcpFactory::createFinalCollectionFromInitial(
- EVENT::LCCollection* mcpVecInitial) {
- // coll must exist
- if (mcpVecInitial) {
- // get num Mcp in coll
- int numInitMcp = mcpVecInitial->getNumberOfElements();
-
- // coll must have members
- if (numInitMcp > 0) {
- // loop over input Mcp coll
- for (int i = 0; i < numInitMcp; i++) {
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "proc initial MCP: " << i << LOG::done;
-#endif
-
- // get Mcp
- MCParticleImpl* mcp =
- static_cast<MCParticleImpl*>(mcpVecInitial->getElementAt(
- i));
-
- // create only from Mcp primaries, e.g. no parents
- if (LcioMcpUtil::isPrimary(mcp)) {
-#ifdef SLIC_LOG
- log() << LOG::debug << "isPrimary" << LOG::done;
-#endif
-
- fromInitialRecurse(mcp);
- }
- }
- } else {
- G4Exception("", "", FatalException,
- "Initial McpVec has no members.");
- }
- } else {
- G4Exception("", "", FatalException, "Initial McpVec ptr is null.");
- }
-}
-
-IMPL::MCParticleImpl* LcioMcpFactory::fromInitialRecurse(EVENT::MCParticle* mcpInit) {
-#ifdef SLIC_LOG
- log() << LOG::debug << "createMcpFromInitialRecurse: " << m_manager->getParticleIndex(m_manager->getInitialMcpCollection(), mcpInit) << LOG::done;
-#endif
-
- MCParticleImpl* mcp = 0;
-
- // find associated primary
- G4PrimaryParticle* g4primary = m_manager->getMaps()->findPrimary(
- mcpInit);
-
- // no associated primary with this initial mcp?
- if (!g4primary) {
-#ifdef SLIC_LOG
- log() << LOG::debug << "initialOnly" << LOG::done;
-#endif
-
- // create from initial Mcp only
- mcp = fromInitialOnly(mcpInit);
- }
- // create from primary
- else {
-#ifdef SLIC_LOG
- log() << LOG::debug << "fromPrimary" << LOG::done;
-#endif
-
- // create it from primary
- mcp = fromPrimary(g4primary, mcpInit);
- }
-
- // return the output MCP
- return mcp;
-}
-
-IMPL::MCParticleImpl* LcioMcpFactory::fromInitialOnly(
- EVENT::MCParticle* mcpInit) {
- log() << LOG::debug << "createMcpFromInitialOnly: "
- << m_manager->getParticleIndex(
- m_manager->getInitialMcpCollection(), mcpInit) << LOG::done;
-
- // copy shallow the initial
- MCParticleImpl* mcp = shallowCopy(mcpInit);
-
- // should add to Mcp map here?
- m_manager->getMaps()->addInitialParticleToFinal(mcpInit, mcp);
-
- // given initial, create and add daughters
- addDaughtersFromInitial(mcp, mcpInit);
-
- // add to final coll
- m_finalColl->addElement(mcp);
-
- return mcp;
-}
-
-IMPL::MCParticleImpl* LcioMcpFactory::fromPrimary(
- G4PrimaryParticle* primary, EVENT::MCParticle* mcpInit) {
-#ifdef SLIC_LOG
- log() << LOG::done << "createMcpFromPrimary - MCP idx <" << m_manager->getParticleIndex(m_manager->getInitialMcpCollection(), mcpInit) << ">" << LOG::done;
-#endif
-
- int trkID = primary->GetTrackID();
-
-#ifdef SLIC_LOG
- if (trkID <= 0)
- {
- log() << "WARNING: trkID <" << trkID << "> looks bad" << LOG::done;
- }
- log() << LOG::done << "primary TID <" << trkID << ">" << LOG::done;
-#endif
-
- MCParticleImpl* mcp = 0;
- Trajectory* trj = TrajectoryManager::instance()->findTrajectory(trkID);
-
- // No trajectory was created.
- if (!trj) {
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "initialAndPrimary" << LOG::done;
-#endif
-
- // Recursively create particles from primary and initial MCParticles.
- mcp = fromInitialAndPrimary(primary, mcpInit);
-
- }
- // Found a trajectory.
- else {
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "initialAndTrajectory" << LOG::done;
-#endif
-
- // Recursively create MCParticles from a trajectory and the initial MCParticle.
- mcp = fromInitialAndTrajectory(trj, mcpInit);
-
- }
-
- // Set generator status from initial particle.
- mcp->setGeneratorStatus(mcpInit->getGeneratorStatus());
-
- // Primary particles do not originate in the simulation.
- mcp->setCreatedInSimulation(false);
-
- // Link the input MCP to an output MCP.
- m_manager->getMaps()->addInitialParticleToFinal(mcpInit, mcp);
-
-//#ifdef SLIC_LOG
-// log() << LOG::debug << "createMcpFromPrimary() - adding trackId <" << trj->GetTrackID() << "> to MCP <" << mcp << "> link." << LOG::done;
-//#endif
-
-// Link the track ID to the output particle.
- if (trkID > 0) {
- m_manager->getMaps()->addTrackToMCParticle(trkID, mcp);
- } else {
-#ifdef SLIC_LOG
- log() << LOG::debug << "WARNING: not adding trkID <" << trkID << "> to map!" << LOG::done;
-#endif
- }
-
-//#ifdef SLIC_LOG
-// log() << LOG::done;
-//#endif
-
- return mcp;
-}
-
-IMPL::MCParticleImpl* LcioMcpFactory::fromInitialAndPrimary(
- G4PrimaryParticle* primary, EVENT::MCParticle* mcpInit) {
- log() << LOG::done << "createMcpFromInitialAndPrimary: "
- << m_manager->getParticleIndex(
- m_manager->getInitialMcpCollection(), mcpInit) << LOG::done;
-
- /*
- * Create new Mcp from primary particle,
- * shallow copy only (no recursion).
- */
- MCParticleImpl* mcpNew = fromPrimaryShallowCopy(primary);
-
- /* Get vertex from initial particle, because primary does not have it -- was immediately decayed */
- double vtx[3] = { mcpInit->getVertex()[0], mcpInit->getVertex()[1],
- mcpInit->getVertex()[2] };
- mcpNew->setVertex(vtx);
-
- /* setup primary's daughters */
- daughtersFromPrimary(primary, /* this primary */
- mcpInit, /* associated initial Mcp particle */
- mcpNew /* parent Mcp particle */
- );
-
-#ifdef SLIC_DEBUG
- if (mcpNew->getDaughters().size() == 0)
- {
-#ifdef SLIC_LOG
- log() << LOG::debug << "No Mcp daughters added." << LOG::done;
-#endif
- }
-#endif
-
- // Copy spin and color flow information if LCIO version is compatible.
-#if LCIO_VERSION_GE(1, 60)
- copySpinAndColorFlow(mcpInit, mcpNew);
-#endif
-
- // done with particle; add to coll
- m_finalColl->addElement(mcpNew);
-
- return mcpNew;
-}
-
-IMPL::MCParticleImpl* LcioMcpFactory::fromInitialAndTrajectory(
- Trajectory* trj, EVENT::MCParticle* mcpInit) {
-#ifdef SLIC_LOG
- log() << LOG::done << "createMcpFromInitialAndTrajectory: " << m_manager->getParticleIndex(m_manager->getInitialMcpCollection(), mcpInit) << LOG::done;
-#endif
-
- // shallow copy create with basic info and no daughters
- IMPL::MCParticleImpl* mcp = fromTrajectory(trj);
-
- // set generator status and created in sim flag
- LcioMcpStatusSetter::setGeneratorStatus(mcpInit, mcp);
-
- // loop over all trajectories to find trajectory daughter
- int numTrj = m_currentTrajectoryContainer->entries();
-
-#ifdef SLIC_LOG
- log() << LOG::done << "nTrajectoryDau: " << numTrj << LOG::done;
-#endif
-
-//#ifdef SLIC_LOG
-// if ( thisTrj->GetParentID() == -1 )
-// {
-// log() << LOG::always << "The trajectory with trackID <" << trj->GetTrackID() << "> has a parent with trackID of -1 !" << LOG::done;
-// }
-//#endif
-
- for (int j = 0; j < numTrj; j++) {
-
- Trajectory* thisTrj =
- static_cast<Trajectory*>((*m_currentTrajectoryContainer)[j]);
-
- // found daughter
- if (thisTrj->GetParentID() == trj->GetTrackID()) {
-#ifdef SLIC_LOG
- log() << LOG::debug << "Trajectory with trackID <" << trj->GetTrackID() << "> has parent with trackID <" << thisTrj->GetParentID() << ">" << LOG::done;
-#endif
-
- // find matching primary Mcp input particle
- MCParticle* mcpPrim =
- m_manager->getMaps()->findInitialParticle(
- thisTrj);
-
- // create MCP output particle from this input daughter
- MCParticleImpl* mcpDauNew = fromInitialAndTrajectory(
- thisTrj, mcpPrim);
-
- // set parent link of new daughter
- mcpDauNew->addParent(mcp);
- }
- // DEBUG
- //else
- //{
- // log() << LOG::done << "BAD NEWS! Could not find Trajectory <" << j << "> with TID <" << trj->GetTrackID() << ">" << LOG::done;
- //}
- // END DEBUG
- }
-
- // add to collection
- m_finalColl->addElement(mcp);
-
- return mcp;
-}
-
-IMPL::MCParticleImpl* LcioMcpFactory::shallowCopy(
- EVENT::MCParticle* mcp) {
- IMPL::MCParticleImpl* mcpNew = new IMPL::MCParticleImpl();
-
- mcpNew->setPDG(mcp->getPDG());
-
- float p[3] = { mcp->getMomentum()[0], mcp->getMomentum()[1],
- mcp->getMomentum()[2] };
-
- mcpNew->setMomentum(p);
- mcpNew->setMass(mcp->getMass());
-
- double vtx[3] = { mcp->getVertex()[0], mcp->getVertex()[1],
- mcp->getVertex()[2] };
- mcpNew->setVertex(vtx);
-
- mcpNew->setGeneratorStatus(mcp->getGeneratorStatus());
- mcpNew->setSimulatorStatus(0);
- mcpNew->setCharge(mcp->getCharge());
- mcpNew->setTime(mcp->getTime());
-
- copySpinAndColorFlow(mcp, mcpNew);
-
- return mcpNew;
-}
-
-// FIXME: Need to pass along colorflow and spin in this method???
-IMPL::MCParticleImpl* LcioMcpFactory::fromPrimaryShallowCopy(
- G4PrimaryParticle* primary) {
- // new MCP
- MCParticleImpl* mcp = new MCParticleImpl();
-
- // PDG
- mcp->setPDG(primary->GetPDGcode());
-
- // momentum
- G4ThreeVector pVec = primary->GetMomentum() / GeV;
- float p[3] = { pVec.x(), pVec.y(), pVec.z() };
- mcp->setMomentum(p);
-
- // mass
- mcp->setMass(primary->GetMass() / GeV);
-
-#ifdef SLIC_LOG
- log() << "mass from primary: " << mcp->getMass() << LOG::endl;
-#endif
-
- // FIXME: Wrong if G4PrimaryParticle actually has valid charge. (Can it happen?)
- mcp->setCharge(LcioMcpManager::m_NAN);
-
- return mcp;
-}
-
-// create LCIO MCParticles from daughters, setting proper linkage
-void LcioMcpFactory::daughtersFromPrimary(G4PrimaryParticle* primary, EVENT::MCParticle* initial, IMPL::MCParticleImpl* parentParticle) {
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "createDaughtersFromPrimary" << LOG::done;
- log() << LOG::debug << "primary: " << primary << LOG::done;
- log() << LOG::debug << "trkID: " << primary->GetTrackID() << LOG::done;
- log() << LOG::debug << "mcpInit: " << initial << LOG::done;
- log() << LOG::debug << "parentParticle: " << parentParticle << LOG::done;
-#endif
-
- // loop over primary daughters
- G4PrimaryParticle* primaryDaughter = primary->GetDaughter();
- while (primaryDaughter) {
-
- // find input Mcp for the primary dau
- MCParticle* mcpDau = m_manager->getMaps()->findDaughterParticle(initial, primaryDaughter);
-
- // must find it or die
- if (0 == mcpDau) {
- G4Exception("", "", FatalException, "Mcp daughter was not found.");
- }
-
- // create dau Mcp
- MCParticleImpl* mcpDauNew = fromPrimary(primaryDaughter, mcpDau);
-
- // make mcpPar the parent
- mcpDauNew->addParent(parentParticle);
-
- // get next primary dau
- primaryDaughter = primaryDaughter->GetNext();
- }
-}
-
-void LcioMcpFactory::addDaughtersFromInitial(IMPL::MCParticleImpl* mcpNew,
- EVENT::MCParticle* mcpInit) {
-#ifdef SLIC_LOG
- log() << LOG::debug << "addMcpDaughtersFromInitial: " << m_manager->getParticleIndex(m_manager->getInitialMcpCollection(), mcpInit) << LOG::done;
-#endif
-
- int numDau = mcpInit->getDaughters().size();
-
- for (int i = 0; i < numDau; i++) {
- // child of input particle
- MCParticle* mcpChildInit = mcpInit->getDaughters()[i];
-
- // use a search function to avoid insert when null!
- MCParticleImpl* mcpChildFinal =
- m_manager->getMaps()->findFinalParticle(
- mcpChildInit);
-
- // if not added already
- if (mcpChildFinal == 0) {
- // create daughter
- mcpChildFinal = fromInitialRecurse(mcpChildInit);
- }
-
- // add parent regardless of above; maybe new MCP or already created
- mcpChildFinal->addParent(mcpNew);
- }
-}
-
-void LcioMcpFactory::fillEndPointEnergy(LCCollectionVec* mcpColl) {
- /* Add endpoint energies to LCGenericObject collection. */
- LcioMcpMaps* maps = LcioMcpManager::instance()->getMaps();
- LCCollectionVec* epColl = new LCCollectionVec(LCIO::LCGENERICOBJECT);
- for (LCCollectionVec::iterator it = mcpColl->begin(); it != mcpColl->end(); it++) {
-
- MCParticle* mcp = static_cast<MCParticle*>(*it);
- G4int trkID = maps->findTrack(mcp);
- double epE = -1.0;
- if (trkID != -1) {
- Trajectory* trj = TrajectoryManager::instance()->findTrajectory(
- trkID);
-
- if (0 != trj) {
- epE = trj->getEndPointEnergy();
- }
- }
-
- LCGenericObjectImpl* obj = new LCGenericObjectImpl();
- obj->setFloatVal(0, epE);
- epColl->push_back(obj);
- }
- LcioManager::instance()->getCurrentLCEvent()->addCollection(epColl, "MCParticleEndPointEnergy");
-}
-
-void LcioMcpFactory::copySpinAndColorFlow(MCParticle* mcp1, MCParticle* mcp2) {
- // Set spin.
- const float* spin1 = mcp1->getSpin();
- float spin2[3];
- spin2[0] = spin1[0];
- spin2[1] = spin1[1];
- spin2[2] = spin1[2];
- ((IMPL::MCParticleImpl*) mcp2)->setSpin(spin2);
-
- // Set color flow.
- const int* cf1 = mcp1->getColorFlow();
- int cf2[2];
- cf2[0] = cf1[0];
- cf2[1] = cf1[1];
- ((IMPL::MCParticleImpl*) mcp2)->setColorFlow(cf2);
-}
-}
projects/slic/trunk/src
--- projects/slic/trunk/src/LcioMcpFilter.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LcioMcpFilter.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,173 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/LcioMcpFilter.cc,v 1.7 2013-11-06 00:23:37 jeremy Exp $
-#include "LcioMcpFilter.hh"
-
-// slic
-#include "LcioMcpPrinter.hh"
-
-// lcio
-#include "IMPL/LCCollectionVec.h"
-#include "EVENT/LCIO.h"
-
-// stl
-#include <cassert>
-
-using EVENT::LCIO;
-
-namespace slic {
-
-LcioMcpFilter::LcioMcpFilter() :
- Module("LcioMcpFilter") {
- ;
-}
-
-LcioMcpFilter::~LcioMcpFilter() {
- ;
-}
-
-EVENT::LCCollection* LcioMcpFilter::filterSimParticles(EVENT::LCCollection* mcpColl) {
- assert(mcpColl);
-
- // new coll w/ gen parts only
- EVENT::LCCollection* genMcpColl = new IMPL::LCCollectionVec(LCIO::MCPARTICLE);
-
- // iterate over input collection
- int mcpSize = mcpColl->getNumberOfElements();
- for (int i = 0; i < mcpSize; i++) {
-
-#ifdef SLIC_DEBUG
- log() << LOG::debug << "handling mcp <" << i << ">" << LOG::done;
-#endif
-
- // next input Mcp
- EVENT::MCParticle* mcp = static_cast<EVENT::MCParticle*>(mcpColl->getElementAt(i));
-
- // primary e.g. no parents?
- if (0 == mcp->getParents().size()) {
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "handling primary Mcp" << LOG::done;
-#endif
-
- // make new Mcp (function checks if gen part)
- makeParticle(mcp, genMcpColl);
-
- // add to coll
- //if ( simMcp ) {
- // genMcpColl->addElement( simMcp );
- //}
- }
- }
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "input coll size <" << mcpSize << ">" << LOG::done;
- log() << LOG::debug << "gen coll size <" << genMcpColl->getNumberOfElements() << ">" << LOG::done;
- log() << LOG::debug << "generator mcp coll..." << LOG::done;
-#endif
-
- //LcioMcpPrinter* prn = new LcioMcpPrinter();
- //prn->printMcpCollection( genMcpColl );
- //delete prn;
-
- m_mcpList.clear();
-
- return genMcpColl;
-}
-
-IMPL::MCParticleImpl* LcioMcpFilter::makeParticle(EVENT::MCParticle* mcp, EVENT::LCCollection* mcpColl) {
- IMPL::MCParticleImpl* mcpGen = 0;
-
- // has non-0 gen status and was not added?
- if (mcp->getGeneratorStatus() && !wasAdded(mcp)) {
-
- //log() << LOG::debug << "gen part" << LOG::done;
-
- // new Mcp part
- mcpGen = new IMPL::MCParticleImpl();
-
- // PDG
- mcpGen->setPDG(mcp->getPDG());
-
- // gen stat
- mcpGen->setGeneratorStatus(mcp->getGeneratorStatus());
-
- // mass
- mcpGen->setMass(mcp->getMass());
-
- // sim stat
- mcpGen->setSimulatorStatus(mcp->getSimulatorStatus());
-
- // momentum
- float p[3] = { mcp->getMomentum()[0], mcp->getMomentum()[1], mcp->getMomentum()[2] };
- mcpGen->setMomentum(p);
-
- // vtx
- double vtx[3] = { mcp->getVertex()[0], mcp->getVertex()[1], mcp->getVertex()[2] };
- mcpGen->setVertex(vtx);
-
- // endp
- double endp[3] = { mcp->getEndpoint()[0], mcp->getEndpoint()[1], mcp->getEndpoint()[2] };
- mcpGen->setEndpoint(endp);
-
- // time
- mcpGen->setTime(mcp->getTime());
-
- // charge
- mcpGen->setCharge(mcp->getCharge());
-
- // add it
- mcpColl->addElement(mcpGen);
- m_mcpList.push_back(mcp);
-
- // add sim daughters
- addSimDaughters(mcp, mcpGen, mcpColl);
- }
- //else {
- // log() << LOG::debug << "skipping sim part and its daughters" << LOG::done;
- //}
-
- return mcpGen;
-}
-
-void LcioMcpFilter::addSimDaughters(EVENT::MCParticle* mcp, EVENT::MCParticle* parMcp,
- EVENT::LCCollection* mcpColl) {
- int nmcp = mcp->getNumberOfDaughters();
-
- //log() << LOG::debug << "handling <" << nmcp << "> dau parts" << LOG::done;
-
- // has at least 1 dau?
- if (nmcp) {
-
- // iterate over dau parts
- for (int i = 0; i < nmcp; i++) {
- EVENT::MCParticle* dauMcp = mcp->getDaughter(i);
-
- // dau != null? (safety check)
- if (dauMcp) {
-
- // make a new dau (could be null )
- IMPL::MCParticleImpl* newDauMcp = makeParticle(dauMcp, mcpColl);
-
- // made a new dau?
- if (newDauMcp) {
- //log() << LOG::debug << "add gen dau" << LOG::done;
-
- // set the dau's par
- newDauMcp->addParent(parMcp);
- }
- }
- }
- }
-}
-
-bool LcioMcpFilter::wasAdded(EVENT::MCParticle* mcp) {
- bool f = false;
- for (std::vector<EVENT::MCParticle*>::const_iterator iter = m_mcpList.begin(); iter != m_mcpList.end();
- iter++) {
- if ((*iter) == mcp) {
- f = true;
- break;
- }
- }
- return f;
-}
-}
projects/slic/trunk/src
--- projects/slic/trunk/src/LcioMcpManager.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LcioMcpManager.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,264 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/LcioMcpManager.cc,v 1.75 2013-11-06 00:23:37 jeremy Exp $
-#include "LcioMcpManager.hh"
-
-// slic
-#include "LcioMcpMessenger.hh"
-#include "LcioManager.hh"
-#include "LcioMcpFactory.hh"
-#include "LcioMcpPrinter.hh"
-#include "TrajectoryManager.hh"
-#include "PrimaryGeneratorAction.hh"
-#include "EventSourceManager.hh"
-
-// lcdd
-#include "lcdd/util/StringUtil.hh"
-
-// lcio
-#include "IMPL/LCCollectionVec.h"
-
-// geant4
-#include "G4Event.hh"
-#include "G4TrajectoryContainer.hh"
-
-// std
-#include <sstream>
-#include <cmath>
-
-using IMPL::MCParticleImpl;
-using EVENT::LCCollection;
-using EVENT::MCParticle;
-using EVENT::LCIO;
-using IMPL::LCCollectionVec;
-
-namespace slic {
-
-const double LcioMcpManager::m_NAN = std::sqrt((double) -1.0);
-const double LcioMcpManager::DEFAULT_MIN_TRACKING_DISTANCE = 1.0;
-
-LcioMcpManager::~LcioMcpManager() {
-}
-
-LcioMcpManager::LcioMcpManager() :
- Module("LcioMcpManager"), m_minimumTrackingDistance(DEFAULT_MIN_TRACKING_DISTANCE), m_enablePrintFinal(
- false), m_enablePrintInitial(false), m_writeInitialMCParticleCollection(false) {
- // create new Mcp messenger instance
- m_messenger = new LcioMcpMessenger();
-
- // setup mcp printer
- m_printer = new LcioMcpPrinter();
-
- /* Setup the MCParticle factory. */
- m_factory = new LcioMcpFactory(this);
-
- /* Create the MCParticle data maps object. */
- m_maps = new LcioMcpMaps(this);
-
- /* Check NAN != NAN */
- assert( m_NAN != m_NAN);
-}
-
-void LcioMcpManager::endEvent(const G4Event* event) {
-#ifdef SLIC_LOG
- log() << LOG::debug << "LcioMcpManager::endEvent" << LOG::done;
-#endif
- /*
- // DEBUG - dump traj cont
- G4TrajectoryContainer* trajectories = event->GetTrajectoryContainer();
- if ( trajectories )
- {
- G4int ntrajectories = trajectories->entries();
-
- if (ntrajectories == 0)
- log() << LOG::debug << "Trajectories container is empty for this event!" << LOG::done;
-
- Trajectory* trj = 0;
- for ( int i=0; i < ntrajectories; i++ )
- {
- trj = static_cast<Trajectory*> ( (*trajectories)[i] );
- #ifdef SLIC_LOG
- log() << LOG::debug << "trajectoryIndex: " << i << LOG::done;
- log() << LOG::debug << "TID: " << trj->GetTrackID() << LOG::done;
-
- //log() << "Trajectory (Name TID PID time): " << ((int)trj->GetTrackID()) << " " << ((int)trj->GetParentID()) << " " << ((std::string)trj->GetParticleName()) << " " << log::DONE;
- log() << LOG::done;
- #endif
- }
- }
- else
- {
- #ifdef SLIC_LOG
- log() << LOG::debug << "No trajectories container for this event!" << LOG::done;
- #endif
- }
- // END DEBUG
- */
-
- /*
- * Create the Mcp collection. Factory will
- * retrieve current event from G4EventManager.
- */
- m_factory->createFinalMcpCollection(event);
-
- if (m_enablePrintFinal) {
- printMcpCollection("final");
- }
-
-#ifdef SLIC_DEBUG
- // print out Mcp processing maps
- m_maps->printMaps();
-#endif
-}
-
-EVENT::LCCollection* LcioMcpManager::findParticleCollection(EVENT::LCEvent* event) {
- const std::vector<std::string>* collNames = event->getCollectionNames();
-
- EVENT::LCCollection* mcpColl = 0;
- for (std::vector<std::string>::const_iterator iter = collNames->begin(); iter != collNames->end();
- iter++) {
- EVENT::LCCollection* thisColl = event->getCollection(*iter);
- if (thisColl->getTypeName() == LCIO::MCPARTICLE) {
- mcpColl = thisColl;
- break;
- }
- }
- return mcpColl;
-}
-
-EVENT::LCCollection* LcioMcpManager::findParticleCollection(const std::string& collName) {
- return m_mcpColls[collName];
-}
-
-void LcioMcpManager::setInitialMcpCollection(EVENT::LCCollection* mcpColl) {
- if (mcpColl) {
- m_mcpColls["initial"] = mcpColl;
- } else {
- G4Exception("", "", FatalException, "Initial collection cannot be set to null.");
- }
-}
-
-EVENT::LCCollection* LcioMcpManager::createParticleCollection(const std::string& collName) {
-#ifdef SLIC_LOG
- log() << LOG::debug << "creating mcp coll <" << collName << ">" << LOG::done;
-#endif
-
- LCCollection* coll = 0;
- if (findParticleCollection(collName)) {
- G4Exception("", "", FatalException, "Attempting to create duplicate Mcp collection.");
- } else {
- coll = new LCCollectionVec(LCIO::MCPARTICLE);
- registerMcpCollection(collName, coll);
- }
-
- return coll;
-}
-
-void LcioMcpManager::registerMcpCollection(const std::string& collName, EVENT::LCCollection* mcpColl) {
- if (mcpColl->getTypeName() != LCIO::MCPARTICLE) {
- G4Exception("", "", FatalException, "The LCCollection does not contain MCParticle objects.");
- }
-
- if (m_mcpColls[collName]) {
- log() << LOG::warning << "WARNING: Mcp collection " + collName + " is already registered"
- << LOG::done;
- } else {
- m_mcpColls[collName] = mcpColl;
- }
-}
-
-void LcioMcpManager::deleteInitialMcpCollection() {
- LCCollection* mcpVec = getInitialMcpCollection();
-
- if (0 != mcpVec) {
- delete mcpVec;
- }
-}
-
-void LcioMcpManager::reset() {
-
- // Initial MCParticle collection needs deletion if it wasn't added to the event.
- if (!writeInitialMCParticleCollection()) {
- deleteInitialMcpCollection();
- }
-
- // Clear the map of MCParticle collections.
- m_mcpColls.clear();
-
- // Clear the maps of MCParticle runtime data from this event.
- m_maps->clear();
-
- // NOTE: The final MCParticle collection is deleted by the LCEvent destructor.
-}
-
-// get MCP collection generated at EndOfEvent
-EVENT::LCCollection* LcioMcpManager::getFinalMcpCollection() {
- if (!findParticleCollection("final")) {
- createParticleCollection("final");
- }
-
- return m_mcpColls["final"];
-}
-
-EVENT::LCCollection* LcioMcpManager::getInitialMcpCollection() {
- if (!findParticleCollection("initial")) {
- createParticleCollection("initial");
- }
-
- return m_mcpColls["initial"];
-}
-
-void LcioMcpManager::setMinimumTrackingDistance(double minDist) {
- m_minimumTrackingDistance = minDist;
-}
-
-double LcioMcpManager::getMinimumTrackingDistance() {
- return m_minimumTrackingDistance;
-}
-
-LcioMcpMaps* LcioMcpManager::getMaps() {
- return m_maps;
-}
-
-void LcioMcpManager::enablePrintFinal(bool p) {
- m_enablePrintFinal = p;
-}
-
-void LcioMcpManager::enablePrintInitial(bool p) {
- m_enablePrintInitial = p;
-}
-
-void LcioMcpManager::beginEvent(const G4Event*) {
- if (m_enablePrintInitial) {
- printMcpCollection("initial");
- }
-}
-
-void LcioMcpManager::printMcpCollection(const std::string& collName) {
- LCCollection* coll = findParticleCollection(collName);
-
- if (coll != 0) {
- printMcpCollection(collName, coll);
- }
-#ifdef SLIC_LOG
- else {
- log() << LOG::warning << "LcioMcpManager::printMcpCollection - " + collName + " was not found" << LOG::done;
- }
-#endif
-}
-
-void LcioMcpManager::printMcpCollection(const std::string& collName, EVENT::LCCollection* coll) {
- m_printer->printParticleCollection(collName, coll);
-}
-
-int LcioMcpManager::getParticleIndex(EVENT::LCCollection* collection, EVENT::MCParticle* particle) {
- if (collection->getTypeName() != LCIO::MCPARTICLE)
- return -1;
-
- for (int i = 0, n = collection->getNumberOfElements(); i < n; i++) {
- EVENT::MCParticle *nextParticle = (MCParticle*) collection->getElementAt(i);
- if (nextParticle == particle) {
- return i;
- }
- }
- return -1;
-}
-}
projects/slic/trunk/src
--- projects/slic/trunk/src/LcioMcpMaps.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LcioMcpMaps.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,185 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/LcioMcpMaps.cc,v 1.8 2013-11-06 00:23:37 jeremy Exp $
-#include "LcioMcpMaps.hh"
-
-// slic
-#include "LcioMcpManager.hh"
-#include "Module.hh"
-#include "Trajectory.hh"
-
-using EVENT::MCParticle;
-using IMPL::MCParticleImpl;
-
-namespace slic {
-
-LcioMcpMaps::LcioMcpMaps(LcioMcpManager* manager) :
- Module("LcioMcpMaps") {
- m_manager = manager;
-}
-
-LcioMcpMaps::~LcioMcpMaps() {
-}
-
-void LcioMcpMaps::addTrackToMCParticle(G4int trackID, MCParticleImpl* particle) {
- assert(particle);
- assert(trackID > 0);
-
- m_trackToMcp[trackID] = particle;
-}
-
-MCParticleImpl* LcioMcpMaps::findParticle(G4int trkID) const {
- MCParticleImpl* mcp = 0;
- for (TrackToMcpMap::const_iterator iter = m_trackToMcp.begin(); iter != m_trackToMcp.end(); iter++) {
- if (iter->first == trkID) {
- mcp = iter->second;
- break;
- }
- }
- return mcp;
-}
-
-/* Find the trackID for the final MCParticle. */
-G4int LcioMcpMaps::findTrack(MCParticle* mcp) const {
- G4int trkID = -1;
- for (TrackToMcpMap::const_iterator iter = m_trackToMcp.begin(); iter != m_trackToMcp.end(); iter++) {
- if (iter->second == mcp) {
- trkID = iter->first;
- break;
- }
- }
- return trkID;
-}
-
-MCParticle* LcioMcpMaps::findInitialParticle(Trajectory* trj) const {
- MCParticle* mcpPrim = 0;
-
- G4int trjTrkID = trj->GetTrackID();
-
- /* Get initial collection. */
- EVENT::LCCollection* initMcpColl = m_manager->getInitialMcpCollection();
-
- // loop over input collection
- int numPrePart = initMcpColl->getNumberOfElements();
- for (int k = 0; k < numPrePart; k++) {
-
- // current MCP
- MCParticle* mcpPre = static_cast<MCParticle*>(initMcpColl->getElementAt(k));
-
- // do we have matching primary for this Mcp?
- G4PrimaryParticle* primaryMatch = findPrimary(mcpPre);
- if (primaryMatch) {
-
- // found matching primary to this input trajectory
- if (primaryMatch->GetTrackID() == trjTrkID) {
- mcpPrim = mcpPre;
- break;
- }
-
- }
- }
-
- return mcpPrim;
-}
-
-MCParticle* LcioMcpMaps::findDaughterParticle(MCParticle* mcpInit, G4PrimaryParticle* primDau) const {
- // loop over Mcp input particle's daughters
- MCParticle* mcpDau = 0;
- int numMcpDau = mcpInit->getDaughters().size();
- for (int i = 0; i < numMcpDau; i++) {
- // does ith Mcp dau match this primary?
- if (findPrimary(mcpInit->getDaughters()[i]) == primDau) {
- // found it
- mcpDau = mcpInit->getDaughters()[i];
- break;
- }
- }
-
- // return dau
- return mcpDau;
-}
-
-MCParticleImpl* LcioMcpMaps::findFinalParticle(MCParticle* mcpInit) const {
- MCParticleImpl* mcpFinal = 0;
- for (InitMcpToMcpMap::const_iterator iter = m_initMcpToMcp.begin(); iter != m_initMcpToMcp.end();
- iter++) {
- if (iter->first == mcpInit) {
- mcpFinal = iter->second;
- break;
- }
- }
-
- return mcpFinal;
-}
-
-G4PrimaryParticle* LcioMcpMaps::findPrimary(MCParticle* mcp) const {
- G4PrimaryParticle *pp = 0;
- for (McpToPrimaryMap::const_iterator iter = m_McpToPrimary.begin(); iter != m_McpToPrimary.end();
- iter++) {
- if (iter->first == mcp) {
- pp = iter->second;
- break;
- }
- }
-
- return pp;
-}
-
-void LcioMcpMaps::addParticleToPrimary(MCParticle* mcp, G4PrimaryParticle* primary) {
- m_McpToPrimary[mcp] = primary;
-}
-
-void LcioMcpMaps::addInitialParticleToFinal(MCParticle* initial, MCParticleImpl* final) {
- m_initMcpToMcp[initial] = final;
-}
-
-void LcioMcpMaps::clear() {
- m_trackToMcp.clear();
- m_McpToPrimary.clear();
- m_initMcpToMcp.clear();
-}
-
-void LcioMcpMaps::printMaps() {
- log() << LOG::debug << "Printing maps..." << LOG::done;
-
- printTracksToParticles();
- printParticlesToPrimaries();
- printInitialParticleToFinal();
- // TrajectoryManager::instance()->printTrackIDToTrajectoryMap();
-}
-
-void LcioMcpMaps::printTracksToParticles() {
- log() << LOG::debug << "Printing track to Mcp map..." << LOG::done;
-
- log() << LOG::debug << "*** TrackToMcpMap ***" << LOG::done;
- log() << LOG::debug << "trkID | McpPtr" << LOG::done;
- for (TrackToMcpMap::iterator iter = m_trackToMcp.begin(); iter != m_trackToMcp.end(); iter++) {
- log() << LOG::debug << iter->first << " " << iter->second << LOG::done;
- }
- log() << LOG::debug << LOG::done;
-}
-
-void LcioMcpMaps::printParticlesToPrimaries() {
- log() << LOG::debug << "Printing Mcp to primary map..." << LOG::done;
-
- log() << LOG::debug << LOG::done;
- log() << LOG::debug << "*** McpToPrimaryMap ***" << LOG::done;
- log() << LOG::debug << "McpPtr | PrimaryPtr" << LOG::done;
-
- for (McpToPrimaryMap::iterator iter = m_McpToPrimary.begin(); iter != m_McpToPrimary.end(); iter++) {
- log() << LOG::debug << iter->first << " " << iter->second << LOG::done;
- }
- log() << LOG::debug << LOG::done;
-}
-
-void LcioMcpMaps::printInitialParticleToFinal() {
- log() << LOG::debug << "Printing InitMcp to Mcp map..." << LOG::done;
-
- log() << LOG::debug << LOG::done;
- log() << LOG::debug << "*** McpToMcpMap ***" << LOG::done;
- log() << LOG::debug << "McpPtrInit | McpPtrFinal" << LOG::done;
-
- for (InitMcpToMcpMap::iterator iter = m_initMcpToMcp.begin(); iter != m_initMcpToMcp.end(); iter++) {
- log() << LOG::debug << iter->first << " " << iter->second << LOG::done;
- }
- log() << LOG::debug << LOG::done;
-}
-}
projects/slic/trunk/src
--- projects/slic/trunk/src/LcioMcpMessenger.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LcioMcpMessenger.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,83 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/LcioMcpMessenger.cc,v 1.14 2013-06-26 01:57:02 jeremy Exp $
-
-// SLIC
-#include "LcioMcpMessenger.hh"
-#include "LcioManager.hh"
-#include "LcioMcpManager.hh"
-#include "TrajectoryManager.hh"
-
-// LCDD
-#include "lcdd/util/StringUtil.hh"
-
-// Geant4
-#include "G4UIparameter.hh"
-
-namespace slic {
-
-LcioMcpMessenger::LcioMcpMessenger() {
- defineCommands();
-}
-
-LcioMcpMessenger::~LcioMcpMessenger() {
-}
-
-void LcioMcpMessenger::SetNewValue(G4UIcommand* cmd, G4String newVals) {
- std::istringstream is((const char*) newVals);
- std::string s;
- is >> s;
-
- LcioMcpManager* mgr = LcioMcpManager::instance();
-
- // set minimum tracking distance
- if (cmd == m_minTrackingDistanceCmd) {
- mgr->setMinimumTrackingDistance(m_minTrackingDistanceCmd->GetNewDoubleValue(newVals));
- } else if (cmd == m_printFinalCollection) {
- mgr->enablePrintFinal(m_printFinalCollection->GetNewBoolValue(newVals.c_str()));
- } else if (cmd == m_printInitialCollection) {
- mgr->enablePrintInitial(m_printInitialCollection->GetNewBoolValue(newVals.c_str()));
- }
- // Enable writing of initial MCParticle collection to event.
- else if (cmd == m_writeInitialCollection) {
- mgr->setWriteInitialMCParticleCollection(m_writeInitialCollection->GetNewBoolValue(newVals.c_str()));
- } else {
- G4Exception("LcioMessenger::setNewValue()", "", JustWarning, "Unknown command");
- }
-
- return;
-}
-
-void LcioMcpMessenger::defineCommands() {
- // Mcp dir
- m_mcpDir = new G4UIdirectory("/mcp/");
- m_mcpDir->SetGuidance("LCIO Monte Carlo Particle commands. [SLIC]");
-
- // minimumTrackingDistance
- m_minTrackingDistanceCmd = new G4UIcmdWithADoubleAndUnit("/mcp/minimumTrackingDistance", this);
- m_minTrackingDistanceCmd->SetGuidance(
- "Set the minimum travel distance for an input particle to be tracked in Geant4.");
- m_minTrackingDistanceCmd->SetParameterName("minimumTrackingDistance", false, true);
- m_minTrackingDistanceCmd->SetDefaultUnit("mm");
- m_minTrackingDistanceCmd->SetDefaultValue(0.1);
- m_minTrackingDistanceCmd->SetRange("minimumTrackingDistance > 0.");
-
- // Enable printing of MCParticle collections during the event.
- m_printInitialCollection = new G4UIcmdWithABool("/mcp/printInitialCollection", this);
- m_printInitialCollection->SetGuidance(
- "Enable printing of the initial MCParticle collection during event processing.");
- m_printInitialCollection->SetParameterName("enable", true);
- m_printInitialCollection->SetDefaultValue(true);
-
- // Enable printing of MCParticle collections during the event.
- m_printFinalCollection = new G4UIcmdWithABool("/mcp/printFinalCollection", this);
- m_printFinalCollection->SetGuidance(
- "Enable printing of the final MCParticle collection during event processing.");
- m_printFinalCollection->SetParameterName("enable", true);
- m_printFinalCollection->SetDefaultValue(true);
-
- m_writeInitialCollection = new G4UIcmdWithABool("/mcp/writeInitialCollection", this);
- m_writeInitialCollection->SetGuidance(
- "Write the initial MCParticle collection from StdHep conversion to the LCEvent.");
- m_writeInitialCollection->SetParameterName("enable", true);
- m_writeInitialCollection->SetDefaultValue(false);
-}
-}
projects/slic/trunk/src
--- projects/slic/trunk/src/LcioMcpPrinter.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LcioMcpPrinter.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,338 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/LcioMcpPrinter.cc,v 1.22 2013-11-06 00:23:37 jeremy Exp $
-
-// SLIC
-#include "LcioMcpPrinter.hh"
-
-// LCDD
-#include "lcdd/util/StringUtil.hh"
-
-// LCIO
-#include "IMPL/LCCollectionVec.h"
-
-using EVENT::LCCollection;
-using IMPL::LCCollectionVec;
-using EVENT::MCParticle;
-using EVENT::MCParticleVec;
-
-using namespace std;
-using std::ios;
-
-namespace slic {
-
-const int LcioMcpPrinter::idx_width = 4;
-const int LcioMcpPrinter::ptr_width = 10;
-const int LcioMcpPrinter::pdg_width = 5;
-const int LcioMcpPrinter::par_width = 15;
-const int LcioMcpPrinter::dau_width = 15;
-const int LcioMcpPrinter::mom_width = 20;
-const int LcioMcpPrinter::vtx_width = 20;
-const int LcioMcpPrinter::end_width = 20;
-const int LcioMcpPrinter::time_width = 6;
-const int LcioMcpPrinter::mass_width = 6;
-const int LcioMcpPrinter::chrg_width = 4;
-const int LcioMcpPrinter::e_width = 6;
-const int LcioMcpPrinter::stat_width = 5;
-const int LcioMcpPrinter::flags_width = 16;
-
-LcioMcpPrinter::LcioMcpPrinter() :
- Module("LcioMcpPrinter"), m_coll(0), m_collName("") {
-}
-
-LcioMcpPrinter::~LcioMcpPrinter() {
-}
-
-void LcioMcpPrinter::printParticleCollection(const std::string& collName, EVENT::LCCollection* coll) {
- m_coll = coll;
- m_collName = collName;
- printParticleCollection();
-}
-
-void LcioMcpPrinter::printParticleCollection() {
- head();
-
- if (m_coll != 0) {
-
- // LCCollectionVec* mcpCollVec = static_cast<LCCollectionVec*> ( m_coll );
-
- if (m_coll->getNumberOfElements() > 0) {
-
- log() << LOG::debug << "Collection Name: " << m_collName << LOG::done;
- log() << LOG::debug << "# of Particles: " << m_coll->getNumberOfElements() << LOG::done;
-
- colLabels();
-
- LCCollectionVec* vec = static_cast<LCCollectionVec*>(m_coll);
- for (LCCollectionVec::iterator iter = vec->begin(); iter != vec->end(); iter++) {
- printParticle(static_cast<EVENT::MCParticle*>(*iter));
- }
-
-#ifdef SLIC_LOG
- flagLabels();
-#endif
-
- } else {
- log() << LOG::debug << "COLLECTION IS EMPTY" << LOG::done;
- }
- } else {
- log() << LOG::debug << "NO CURRENT MCPARTICLE COLLECTION" << LOG::done;
- }
-
- foot();
-}
-
-void LcioMcpPrinter::flagLabels() {
- line();
-
- log() << LOG::debug << "FLAGS: createdInSimulation Backscatter vertexIsNotEndpointOfParent "
- << "DecayedInTracker DecayedInCalorimeter LeftDetector Stopped";
-
- line();
-}
-
-void LcioMcpPrinter::printParticle(EVENT::MCParticle* mcp) {
- // index
- log().getOutputStream().width(idx_width);
- log().getOutputStream() << left;
- log() << LOG::debug << findParticleIndex(mcp);
-
- sep();
-
- // ptr
- log().getOutputStream().width(ptr_width);
- log().getOutputStream() << left;
- log() << LOG::debug << mcp;
-
- sep();
-
- // PDG
- log().getOutputStream().width(pdg_width);
- log() << LOG::debug << left << mcp->getPDG();
-
- sep();
-
- // parents
- log().getOutputStream().width(par_width);
- log() << LOG::debug << left << makeIndexVectorString(findParentIndices(mcp));
-
- sep();
-
- // daughters
- log().getOutputStream().width(dau_width);
- log() << LOG::debug << left << makeIndexVectorString(findDaughterIndices(mcp));
-
- sep();
-
- // momentum
- log().getOutputStream().width(mom_width);
- log() << LOG::debug << left << makeDoubleArray3String(mcp->getMomentum());
-
- sep();
-
- // vertex
- log().getOutputStream().width(vtx_width);
- log() << LOG::debug << left << makeDoubleArray3String(mcp->getVertex());
-
- sep();
-
- // endpoint
- log().getOutputStream().width(end_width);
- log() << LOG::debug << left << makeDoubleArray3String(mcp->getEndpoint());
-
- sep();
-
- // time
- log().getOutputStream().width(time_width);
- log() << LOG::debug << left << mcp->getTime();
-
- sep();
-
- // mass
- log().getOutputStream().width(mass_width);
- log() << LOG::debug << left << mcp->getMass();
-
- sep();
-
- // charge
- log().getOutputStream().width(chrg_width);
- log() << LOG::debug << mcp->getCharge();
-
- sep();
-
- // energy
- log().getOutputStream().width(e_width);
- log() << LOG::debug << mcp->getEnergy();
-
- // gen status
- log() << LOG::debug << mcp->getGeneratorStatus();
-
-#ifdef SLIC_LOG
- sep();
- log() << LOG::debug << left << makeSimStatusCodesString( mcp );
-#endif
-
- line();
- line();
-}
-
-std::string LcioMcpPrinter::makeDoubleArray3String(const double* arr) {
-
- return "(" + StringUtil::toString(arr[0]) + ", " + StringUtil::toString(arr[1]) + ", "
- + StringUtil::toString(arr[2]) + ")";
-}
-
-std::string LcioMcpPrinter::makeFloatArray3String(const float* arr) {
- return "(" + StringUtil::toString(arr[0]) + ", " + StringUtil::toString(arr[1]) + ", "
- + StringUtil::toString(arr[2]) + ")";
-}
-
-std::string LcioMcpPrinter::makeSimStatusCodesString(EVENT::MCParticle* mcp) {
- return StringUtil::toString(mcp->isCreatedInSimulation()) + " "
- + StringUtil::toString(mcp->isBackscatter()) + " "
- + StringUtil::toString(mcp->vertexIsNotEndpointOfParent()) + " "
- + StringUtil::toString(mcp->isDecayedInTracker()) + " "
- + StringUtil::toString(mcp->isDecayedInCalorimeter()) + " "
- + StringUtil::toString(mcp->hasLeftDetector()) + " " + StringUtil::toString(mcp->isStopped());
-}
-
-std::string LcioMcpPrinter::makeIndexVectorString(const IndexVector& idxVec) {
- std::string idx_vec;
- idx_vec = "[ ";
-
- for (IndexVector::const_iterator iter = idxVec.begin(); iter != idxVec.end(); iter++) {
- int this_idx = (*iter);
- idx_vec += StringUtil::toString(this_idx) + " ";
- }
-
- idx_vec += "]";
- return idx_vec;
-}
-
-LcioMcpPrinter::IndexVector LcioMcpPrinter::findDaughterIndices(EVENT::MCParticle* mcp) {
- IndexVector dauIdxs;
- const EVENT::MCParticleVec& dauVec = mcp->getDaughters();
- for (EVENT::MCParticleVec::const_iterator iter = dauVec.begin(); iter != dauVec.end(); iter++) {
- dauIdxs.push_back(findParticleIndex((const_cast<MCParticle*>(*iter))));
- }
-
- return dauIdxs;
-}
-
-LcioMcpPrinter::IndexVector LcioMcpPrinter::findParentIndices(EVENT::MCParticle* mcp) {
- IndexVector parIdxs;
- const EVENT::MCParticleVec& parVec = mcp->getParents();
- for (EVENT::MCParticleVec::const_iterator iter = parVec.begin(); iter != parVec.end(); iter++) {
- parIdxs.push_back(findParticleIndex((const_cast<MCParticle*>(*iter))));
- }
-
- return parIdxs;
-}
-
-int LcioMcpPrinter::findParticleIndex(EVENT::MCParticle* mcp) {
- int idx = 0;
- bool fnd = false;
- LCCollectionVec* mcpCollVec = static_cast<LCCollectionVec*>(m_coll);
- for (LCCollectionVec::iterator iter = mcpCollVec->begin(); iter != mcpCollVec->end(); iter++) {
- // ptr compare
- if ((static_cast<MCParticle*>(*iter)) == mcp) {
- fnd = true;
- break;
- }
- idx++;
- }
-
- if (!fnd) {
- idx = -1;
- }
-
- return idx;
-}
-
-inline void LcioMcpPrinter::sep() {
- log() << LOG::debug << "|";
-}
-
-inline void LcioMcpPrinter::line() {
- log() << LOG::debug << LOG::done;
-}
-
-inline void LcioMcpPrinter::head() {
- log() << LOG::debug << "---------------------MCParticle Collection---------------------";
-
- line();
-}
-
-inline void LcioMcpPrinter::foot() {
- log() << LOG::debug << "---------------------------------------------------------------";
-
- line();
-}
-
-void LcioMcpPrinter::colLabels() {
- log().getOutputStream().width(idx_width);
- log() << LOG::debug << left << "idx";
-
- sep();
-
- log().getOutputStream().width(ptr_width);
- log() << LOG::debug << left << "ptr";
-
- sep();
-
- log().getOutputStream().width(pdg_width);
- log() << LOG::debug << left << "pdg";
-
- sep();
-
- log().getOutputStream().width(par_width);
- log() << LOG::debug << left << "par";
-
- sep();
-
- log().getOutputStream().width(dau_width);
- log() << LOG::debug << left << "dau";
-
- sep();
-
- log().getOutputStream().width(mom_width);
- log() << LOG::debug << left << "mom";
-
- sep();
-
- log().getOutputStream().width(vtx_width);
- log() << LOG::debug << left << "vtx";
-
- sep();
-
- log().getOutputStream().width(end_width);
- log() << LOG::debug << left << "end";
-
- sep();
-
- log().getOutputStream().width(time_width);
- log() << LOG::debug << left << "time";
-
- sep();
-
- log().getOutputStream().width(mass_width);
- log() << LOG::debug << left << "mass";
-
- sep();
-
- log().getOutputStream().width(chrg_width);
- log() << LOG::debug << left << "chg";
-
- sep();
-
- log().getOutputStream().width(e_width);
- log() << LOG::debug << "E";
-
- sep();
-
-#ifdef SLIC_LOG
- log().getOutputStream().width(flags_width);
- log() << LOG::debug << "flags";
-#endif
-
- line();
-}
-}
projects/slic/trunk/src
--- projects/slic/trunk/src/LcioMcpStatusSetter.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LcioMcpStatusSetter.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,99 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/LcioMcpStatusSetter.cc,v 1.7 2013-11-06 00:23:37 jeremy Exp $
-#include "LcioMcpStatusSetter.hh"
-
-// slic
-#include "Trajectory.hh"
-#include "EventSourceManager.hh"
-
-using IMPL::MCParticleImpl;
-using EVENT::MCParticle;
-
-namespace slic {
-IMPL::MCParticleImpl* LcioMcpStatusSetter::setStatusCodesFromTrajectory(Trajectory* trj,
- MCParticleImpl* mcp) {
- // set vertexIsNotEndOfParent flag
- LcioMcpStatusSetter::setVertexIsNotEndpointOfParentFlag(trj, mcp);
-
- // set final status flags
- LcioMcpStatusSetter::setFinalStatusFlags(trj, mcp);
-
- // set backscatter flag
- LcioMcpStatusSetter::setBackscatterFlag(trj, mcp);
-
- // For GPS source, set generator status to FS for trajectories without parents,
- // i.e. the parent track ID is 0.
- if (EventSourceManager::instance()->getCurrentSourceType() == EventSourceManager::eGPS) {
- if (trj->GetParentID() == 0) {
- mcp->setGeneratorStatus(1);
- }
- }
-
- // return the MCP
- return mcp;
-}
-
-void LcioMcpStatusSetter::setBackscatterFlag(Trajectory* trj, IMPL::MCParticleImpl* mcp) {
- if (trj->getBackscatter()) {
- mcp->setBackscatter(true);
- }
-}
-
-void LcioMcpStatusSetter::setVertexIsNotEndpointOfParentFlag(Trajectory* trj, MCParticleImpl* mcp) {
- mcp->setVertexIsNotEndpointOfParent(trj->getVertexIsNotEndpointOfParent());
-}
-
-void LcioMcpStatusSetter::setCreatedInSimulationFlag(Trajectory* trj, IMPL::MCParticleImpl* mcp) {
- if (trj->getCreatedInSimulation()) {
- mcp->setGeneratorStatus(0);
- mcp->setCreatedInSimulation(true);
- } else {
- mcp->setCreatedInSimulation(false);
- }
-}
-
-void LcioMcpStatusSetter::setFinalStatusFlags(Trajectory* trj, IMPL::MCParticleImpl* mcp) {
- // set final status flags
- Trajectory::EFinalStatus trjFinalStatus = trj->getFinalStatus();
- switch (trjFinalStatus) {
-
- case Trajectory::eDecayedInTracker:
- mcp->setDecayedInTracker(true);
- break;
-
- case Trajectory::eLeftDetector:
- mcp->setHasLeftDetector(true);
- break;
-
- case Trajectory::eStopped:
- mcp->setStopped(true);
- break;
-
- case Trajectory::eDecayedInCalorimeter:
- mcp->setDecayedInCalorimeter(true);
- break;
-
- case Trajectory::eUnset:
- break;
-
- default:
- break;
- }
-}
-
-void LcioMcpStatusSetter::setGeneratorStatus(MCParticle* mcpInit, MCParticleImpl* mcp) {
- // set generator status if initial exists
- if (mcpInit) {
-
- // if init particle, this was not created in sim
- mcp->setCreatedInSimulation(false);
-
- // copy gen status from initial
- mcp->setGeneratorStatus(mcpInit->getGeneratorStatus());
- }
- // no initial -> created in sim
- else {
- mcp->setCreatedInSimulation(true);
- mcp->setGeneratorStatus(0);
- }
-}
-}
projects/slic/trunk/src
--- projects/slic/trunk/src/LcioMcpUtil.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LcioMcpUtil.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,31 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/LcioMcpUtil.cc,v 1.4 2013-11-06 00:23:37 jeremy Exp $
-#include "LcioMcpUtil.hh"
-
-namespace slic {
-// Mcp with no parents
-bool LcioMcpUtil::isPrimary(IMPL::MCParticleImpl* mcp) {
- return (mcp->getParents().size() == 0);
-}
-
-G4double LcioMcpUtil::computeParticleDistance(IMPL::MCParticleImpl* mcp1, IMPL::MCParticleImpl* mcp2) {
- return computeDistance(mcp1->getVertex(), mcp2->getVertex());
-}
-
-G4double LcioMcpUtil::computeDistance(const double *v1, const double *v2) {
- // distance = sqrt ( x^2 + y^2 + z^2 )
-
- // x^2
- G4double x = (v1[0] - v2[0]);
- G4double x2 = pow(x, 2);
-
- // y^2
- G4double y = (v1[1] - v2[1]);
- G4double y2 = pow(y, 2);
-
- // z^2
- G4double z = (v1[2] - v2[2]);
- G4double z2 = pow(z, 2);
-
- return sqrt(x2 + y2 + z2);
-}
-}
projects/slic/trunk/src
--- projects/slic/trunk/src/LcioMessenger.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LcioMessenger.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -3,9 +3,7 @@
// SLIC
#include "LcioMessenger.hh"
#include "LcioManager.hh"
-#include "LcioMcpManager.hh"
#include "LcioFileNamer.hh"
-#include "TrajectoryManager.hh"
// LCDD
#include "lcdd/util/StringUtil.hh"
@@ -62,8 +60,8 @@
m_mgr->setAutonameFields(fields);
}
/* dump event */
- else if (cmd == m_dumpEventCmd) {
- m_mgr->enableDumpEvent(m_dumpEventCmd->GetNewBoolValue(newVals.c_str()));
+ else if (cmd == m_printEventCmd) {
+ m_mgr->enableDumpEvent(m_printEventCmd->GetNewBoolValue(newVals.c_str()));
}
// flags
else if (cmd == m_setLongFlagCmd || cmd == m_setPDGFlagCmd) {
@@ -76,7 +74,7 @@
flag_set = StringUtil::toBool(s);
}
- LcioHitsCollectionBuilder* hcb = m_mgr->getHCBuilder();
+ LcioHitsCollectionBuilder* hcb = m_mgr->getHitsCollectionBuilder();
if (hcb) {
@@ -92,6 +90,10 @@
G4Exception("", "", FatalException, "LcioHitsCollectionBuilder is null.");
}
}
+ else if (cmd == _writeCompleteEventCmd) {
+ G4bool writeCompleteEvent = StringUtil::toBool(s);
+ LcioManager::instance()->setWriteCompleteEvent(writeCompleteEvent);
+ }
// bad command; shouldn't happen
else {
G4Exception("LcioMessenger::setNewValue()", "", JustWarning, "Unknown command.");
@@ -154,9 +156,14 @@
m_autonameCmd->SetParameter(p);
/* dumping of event data */
- m_dumpEventCmd = new G4UIcmdWithABool("/lcio/dumpEvent", this);
- m_dumpEventCmd->SetGuidance("Dump information about collections in the event");
- m_dumpEventCmd->SetParameterName("enable", true);
- m_dumpEventCmd->SetDefaultValue(true);
+ m_printEventCmd = new G4UIcmdWithABool("/lcio/printEvent", this);
+ m_printEventCmd->SetGuidance("Dump information about collections in the event");
+ m_printEventCmd->SetParameterName("enable", true);
+ m_printEventCmd->SetDefaultValue(true);
+
+ _writeCompleteEventCmd = new G4UIcmdWithABool("/lcio/writeCompleteEvent", this);
+ _writeCompleteEventCmd->SetGuidance("Write complete event to MCParticle output.");
+ _writeCompleteEventCmd->SetParameterName("enable", true);
+ _writeCompleteEventCmd->SetDefaultValue(true);
}
}
projects/slic/trunk/src
--- projects/slic/trunk/src/LcioPrimaryGenerator.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/LcioPrimaryGenerator.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,261 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/LcioPrimaryGenerator.cc,v 1.42 2013-11-06 00:23:37 jeremy Exp $
-
-// SLIC
-#include "LogStream.hh"
-#include "LcioPrimaryGenerator.hh"
-#include "LcioManager.hh"
-#include "LcioMcpManager.hh"
-#include "LcioMcpUtil.hh"
-
-// LCDD
-#include "lcdd/util/StringUtil.hh"
-
-// Geant4
-#include "G4ParticleTable.hh"
-#include "G4SystemOfUnits.hh"
-
-// STL
-#include <sstream>
-
-using IMPL::MCParticleImpl;
-using EVENT::LCCollection;
-using EVENT::LCIO;
-
-namespace slic {
-
-LcioPrimaryGenerator::LcioPrimaryGenerator(LcioManager* mgr) :
- Module("LcioPrimaryGenerator") {
- m_mgr = mgr;
-
- m_mcpManager = LcioMcpManager::instance();
-}
-
-LcioPrimaryGenerator::~LcioPrimaryGenerator() {
-}
-
-void LcioPrimaryGenerator::generatePrimaryVertexFromMcpCollection(LCCollection* mcpVec, G4Event* anEvent) {
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "********** Generating Event from LCIO MCParticles **********" << LOG::endl << LOG::done;
-#endif
-
- assert( mcpVec);
- assert( mcpVec->getTypeName() == LCIO::MCPARTICLE);
-
- G4int nhep = mcpVec->getNumberOfElements();
-
- if (nhep < 1) {
- return;
- }
-
-//#ifdef SLIC_LOG
-// log() << LOG::debug << "nInitialMCParticles: " << nhep << LOG::done;
-//#endif
-
- for (int i = 0; i < nhep; i++) {
-
- MCParticleImpl* mcp = dynamic_cast<MCParticleImpl*>(mcpVec->getElementAt(i));
-
- //assert( mcp );
- if (mcp == 0)
- G4Exception("", "", FatalException, "Could not find MCParticle at indx.");
-
- // debug info
- G4ParticleTable* tbl = G4ParticleTable::GetParticleTable();
- G4ParticleDefinition* pdef = tbl->FindParticle(mcp->getPDG());
-
- G4String pname = "unknown";
- if (pdef != 0)
- pname = pdef->GetParticleName();
- // end debug info
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "generating particle #" << i << ", pdg: " << mcp->getPDG() << ", name: " << pname << LOG::done;
- log() << LOG::debug << "vertex: ( " << mcp->getVertex()[0] << ", " << mcp->getVertex()[1] << ", " << mcp->getVertex()[2] << " )" << LOG::done;
- if (mcp->getEndpoint() != 0)
- log() << LOG::debug << "endp: ( " << mcp->getEndpoint()[0] << ", " << mcp->getEndpoint()[1] << ", " << mcp->getEndpoint()[2] << " )" << LOG::done;
- log() << LOG::debug << "time: " << mcp->getTime() << LOG::done;
-#endif
-
- G4int numPar = mcp->getParents().size();
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "nparents: " << numPar << LOG::done;
-#endif
-
- /*
- * Handle 4 cases, following Ron Cassell's LCSHEPEvtInterface from LCS package:
- * case 1: predecay
- * case 2: not predecay; "final state"
- * case 3: not predecay; "intermediate" or "documentation" TO BE TRACKED
- * case 4: not predecay; "intermediate" or "documentation" *NOT* TRACKED
- */
- // loop vars
- MCParticleImpl* parMcp = 0;
- G4PrimaryParticle* g4parent = 0;
-
- G4PrimaryParticle* thePrimary = 0;
- G4PrimaryVertex* theVertex = 0;
-
- G4bool isPreDecay = false;
- G4bool createPrimary = false;
- G4bool createVertex = false;
-
- // check if case 1 (PreDecay)
- if (numPar > 0) {
- parMcp = dynamic_cast<MCParticleImpl*>(mcp->getParents()[0]);
- g4parent = m_mcpManager->getMaps()->findPrimary(parMcp);
-
- // Check if the particle has parents, in which case it needs a predecay.
- if (g4parent != 0) {
-
- isPreDecay = true;
-#ifdef SLIC_LOG
- log() << "PREDECAY" << LOG::done;
- log() << LOG::debug << "parent idx: " << LcioMcpManager::instance()->getParticleIndex(mcpVec, parMcp) << LOG::done;
-#endif
- }
- }
-
- // case 1
- if (isPreDecay) {
- createPrimary = true;
- // no vertex
- } else {
- // case 2: final state
- if (mcp->getGeneratorStatus() == 1) {
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "FINAL STATE" << LOG::done;
-#endif
-
- createPrimary = true;
- createVertex = true;
- } else {
- G4double dist = 0;
-
- // check case 3, if traveled > minDist
- if (mcp->getDaughters().size() > 0) {
-
- IMPL::MCParticleImpl* firstDau =
- dynamic_cast<IMPL::MCParticleImpl*>(mcp->getDaughters()[0]);
-
- dist = LcioMcpUtil::computeParticleDistance(mcp, firstDau);
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "mcpDistance=" << dist << LOG::done;
- log() << LOG::debug << "minTrackingDist=" << m_mcpManager->getMinimumTrackingDistance() << LOG::done;
-#endif
-
- // case 3
- if (dist > m_mcpManager->getMinimumTrackingDistance()) {
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "INTERMEDIATE OR DOC TO BE TRACKED" << LOG::done;
-#endif
-
- createPrimary = true;
- createVertex = true;
- }
- // case 4 = no-op
- // *This particle will not be tracked.*
-#ifdef SLIC_LOG
- else {
- log() << LOG::debug << "INTERMEDIATE OR DOC THAT WILL NOT BE TRACKED" << LOG::done;
- }
-#endif
- }
-
- }
- }
-
-#ifdef SLIC_LOG
- log() << "isPreDecay: " << isPreDecay << LOG::done;
- log() << "createPrimary: " << createPrimary << LOG::done;
- log() << "createVertex: " << createVertex << LOG::done;
-#endif
-
- // create a primary
- if (createPrimary) {
-
- thePrimary = createPrimaryParticle(mcp);
-
- // set daughters for PreDecay
- if (isPreDecay) {
-
- assert( g4parent);
-
- //
- // Computation of proper_time from RC
- //
-
- G4ThreeVector parMom = g4parent->GetMomentum();
-
-#ifdef SLIC_LOG
- log() << "parent p: " << parMom << LOG::done;
-#endif
-
- G4double E = sqrt(
- pow(g4parent->GetMass(), 2) + pow(parMom.x(), 2) + pow(parMom.y(), 2)
- + pow(parMom.z(), 2));
-
-#ifdef SLIC_LOG
- log() << "particle time: " << mcp->getTime() << LOG::done;
- log() << "parent time: " << parMcp->getTime() << LOG::done;
- log() << "parent mass: " << g4parent->GetMass() << LOG::done;
-#endif
-
- G4double proper_time = ((mcp->getTime() - parMcp->getTime()) * g4parent->GetMass()) / E;
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "parent E: " << E << LOG::done;
- log() << LOG::debug << "proper time: " << proper_time << LOG::done;
-#endif
-
- g4parent->SetDaughter(thePrimary);
- g4parent->SetProperTime(proper_time);
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "mcp decay time: " << mcp->getTime() - parMcp->getTime() << LOG::done;
-#endif
- }
- }
-
- // create a vertex, add primary and set in event
- if (createVertex) {
- theVertex = createPrimaryVertex(mcp);
- theVertex->SetPrimary(thePrimary);
- anEvent->AddPrimaryVertex(theVertex);
- }
-
- // insert mcp, primary pair into LcioManager's map (could be null)
- if (thePrimary) {
- m_mcpManager->getMaps()->addParticleToPrimary(mcp, thePrimary);
- }
-
-#ifdef SLIC_LOG
- log() << LOG::debug << LOG::endl;
-#endif
- }
-}
-
-G4PrimaryParticle* LcioPrimaryGenerator::createPrimaryParticle(IMPL::MCParticleImpl* mcp) {
- G4PrimaryParticle* primary = new G4PrimaryParticle(mcp->getPDG(), mcp->getMomentum()[0] * GeV,
- mcp->getMomentum()[1] * GeV, mcp->getMomentum()[2] * GeV);
- primary->SetMass(mcp->getMass() * GeV);
-
- return primary;
-}
-
-G4PrimaryVertex* LcioPrimaryGenerator::createPrimaryVertex(IMPL::MCParticleImpl* mcp) {
-#ifdef SLIC_LOG
- log() << "createPrimaryVertexFromMcp" << LOG::done;
- log() << "particle time: " << mcp->getTime() << LOG::done;
-#endif
-
- G4ThreeVector pos = G4ThreeVector(mcp->getVertex()[0], mcp->getVertex()[1], mcp->getVertex()[2]);
-
- return new G4PrimaryVertex(pos, mcp->getTime());
-
-}
-}
projects/slic/trunk/src
--- projects/slic/trunk/src/ParticleGunEventSource.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/ParticleGunEventSource.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -19,7 +19,7 @@
m_generator->GeneratePrimaryVertex(evt);
}
-void ParticleGunEventSource::dumpCurrentEvent() {
+void ParticleGunEventSource::printCurrentEvent() {
log() << LOG::okay << "Dumping ParticleGun Event info..." << LOG::done;
G4ParticleDefinition* pdef = m_generator->GetParticleDefinition();
projects/slic/trunk/src
--- projects/slic/trunk/src/PrimaryGeneratorAction.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/PrimaryGeneratorAction.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -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,14 +16,11 @@
#include "IMPL/MCParticleImpl.h"
#include "EVENT/LCCollection.h"
-// Geant4
-#include "G4Types.hh"
-#include "G4StateManager.hh"
-#include "G4SystemOfUnits.hh"
-
namespace slic {
+
PrimaryGeneratorAction::PrimaryGeneratorAction() :
Module("PrimaryGeneratorAction", false) {
+ _manager = EventSourceManager::instance();
}
PrimaryGeneratorAction::~PrimaryGeneratorAction() {
@@ -31,161 +28,34 @@
void PrimaryGeneratorAction::GeneratePrimaries(G4Event *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;
- }
+ /* Print begin of event message. */
+ printBeginEventMessage(anEvent);
- EventSourceManager* mgr = EventSourceManager::instance();
+ /* Reset the TrackManager's state from previous event. */
+ TrackManager::instance()->reset();
- // Beginning of event message that prints the event number.
- // This should show up as the first line in the log messages for this event.
- printBeginEventMessage(anEvent);
+ /* Reset the MCParticleManager's state from previous event. */
+ MCParticleManager::instance()->reset();
- // State of MCP manager needs to be reset for new event.
- // FIXME: This is also called in LcioManager::reset() but not working for some reason!
- LcioMcpManager::instance()->reset();
+ /* 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 which will read the next input event if applicable. */
+ _manager->beginEvent(anEvent);
- // Generator manager begin event hook.
- mgr->beginEvent(anEvent);
+ /* Generate primary vertices using the event generator's (previously read) current event. */
+ _manager->generateNextEvent(anEvent);
+ }
- // Generate vertices using the manager.
- mgr->GeneratePrimaryVertex(anEvent);
-
- // 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;
}
-
-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/trunk/src
--- projects/slic/trunk/src/RunAction.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/RunAction.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,9 +1,11 @@
// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/RunAction.cc,v 1.23 2013-06-26 01:57:02 jeremy Exp $
// SLIC
+//#include "LcioManager.hh"
+#include "EventSourceManager.hh"
#include "LcioManager.hh"
-#include "EventSourceManager.hh"
-#include "EventDebugger.hh"
+//#include "EventDebugger.hh"
+#include "TrackManager.hh"
// LCDD
#include "lcdd/core/GeometryManager.hh"
@@ -22,6 +24,7 @@
}
void RunAction::BeginOfRunAction(const G4Run *aRun) {
+
// Start the run clock.
startRunTimer();
@@ -44,6 +47,7 @@
}
void RunAction::EndOfRunAction(const G4Run *run) {
+
// LcioManager's action
LcioManager::instance()->endRun(run);
projects/slic/trunk/src
--- projects/slic/trunk/src/RunManager.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/RunManager.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -9,6 +9,7 @@
#include "PrimaryGeneratorAction.hh"
#include "RunAction.hh"
#include "SlicApplication.hh"
+#include "StackingAction.hh"
#include "SteppingAction.hh"
#include "TrackingAction.hh"
@@ -30,11 +31,12 @@
}
void RunManager::initializeUserActions() {
- SetUserAction(new PrimaryGeneratorAction);
- SetUserAction(new RunAction);
- SetUserAction(new EventAction);
- SetUserAction(new TrackingAction);
- SetUserAction(new SteppingAction);
+ SetUserAction(new PrimaryGeneratorAction());
+ SetUserAction(new RunAction());
+ SetUserAction(new EventAction());
+ SetUserAction(new TrackingAction());
+ SetUserAction(new SteppingAction());
+ SetUserAction(new StackingAction());
m_userActionsInitialized = true;
}
projects/slic/trunk/src
--- projects/slic/trunk/src/SlicApplication.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/SlicApplication.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -4,7 +4,6 @@
#include "SlicApplication.hh"
#include "CommandLineProcessor.hh"
#include "CommandQueue.hh"
-#include "EventDebugger.hh"
#include "FieldManager.hh"
#include "SlicApplicationMessenger.hh"
#include "LcioManager.hh"
@@ -12,38 +11,29 @@
#include "PackageInfo.hh"
#include "RunManager.hh"
#include "TimeUtil.hh"
-#include "TrajectoryManager.hh"
#include "HepPDTManager.hh"
#include "VRML2WriterMessenger.hh"
+#include "PhysicsListManager.hh"
// LCDD
#include "lcdd/util/StringUtil.hh"
#include "lcdd/core/LCDDDetectorConstruction.hh"
-// geant4
+// Geant4
#include "G4ApplicationState.hh"
-#include "G4ScoringManager.hh"
#include "G4StateManager.hh"
-// geant4 visualization
#ifdef G4VIS_USE
#include "G4VisExecutive.hh"
#endif
-// geant4 UI
#include "G4UIExecutive.hh"
-#ifdef HAVE_G4VERSION
-#include "G4Version.hh"
-#endif
-
-#include "PhysicsListManager.hh"
-
namespace slic {
SlicApplication::SlicApplication() :
Module("SlicApplication", false), m_appMessenger(0), m_runManager(0), m_mode(SlicApplication::eBatch), m_returnCode(
- 0), m_setRunAbort(false), m_isInitialized(false), m_geant4VersionString(0) {
+ 0), m_isInitialized(false), m_geant4VersionString(0) {
#ifdef SLIC_LOG
log() << LOG::okay << LOG::head << "SLIC is starting." << LOG::done;
#endif
@@ -101,9 +91,6 @@
// Create a new RunManager.
m_runManager = new RunManager();
- // Activate command-based scorer.
- G4ScoringManager::GetScoringManager();
-
// Initialize the physics list manager.
PhysicsListManager::instance();
@@ -177,19 +164,12 @@
return m_mode;
}
-void SlicApplication::setAborting(bool a) {
- if (a) {
- // FIXME: Does this actually do anything to current G4Event?
- G4RunManager::GetRunManager()->AbortRun(a);
+void SlicApplication::setAborting(bool abort) {
+ if (abort) {
+ G4RunManager::GetRunManager()->AbortRun(abort);
}
-
- m_setRunAbort = a;
}
-bool SlicApplication::isAborting() const {
- return m_setRunAbort;
-}
-
RunManager* SlicApplication::getRunManager() {
return m_runManager;
}
@@ -228,16 +208,8 @@
void SlicApplication::printSplashScreen() {
log() << LOG::done;
log() << LOG::always << "*************************************************************" << LOG::done;
- log() << LOG::always << " Application : " << PackageInfo::getNameString() << LOG::done;
- //" ("
- //<< PackageInfo::getAbbrevString() << ")" << LOG::done;
- log() << LOG::always << " Version : " << PackageInfo::getVersionString() << LOG::done;
-
- //log() << LOG::always << " Date : " << PackageInfo::getChangeDateString() << LOG::done;
- //log() << LOG::always << " Authors : " << PackageInfo::getAuthorString() << LOG::done;
- //log() << LOG::always << " Inst : " << PackageInfo::getInstitutionString() << LOG::done;
- //log() << LOG::always << " WWW : " << PackageInfo::getWWW() << LOG::done;
- //log() << LOG::always << " Contact : " << PackageInfo::getEmail() << LOG::done;
+ log() << LOG::always << " Application : " << PackageInfo::getFullName() << " (" << PackageInfo::getShortName() << ")" << LOG::done;
+ log() << LOG::always << " Version : " << PackageInfo::getVersion() << LOG::done;
log() << LOG::always << "*************************************************************" << LOG::done;
}
@@ -245,31 +217,4 @@
CommandLineProcessor::instance()->printUsage();
}
-std::string SlicApplication::getGeant4VersionString() {
- if (m_geant4VersionString == 0) {
- std::string g4ver;
-#ifdef HAVE_G4VERSION
- g4ver = StringUtil::toString(G4VERSION_NUMBER);
-#else
- // Massage the string returned by G4RunManager into format
- // returned by G4Version in recent Geant4 versions.
- g4ver = G4RunManager::GetRunManager()->GetVersionString();
- std::string::size_type start = g4ver.find("geant4-", 0) + 7;
- std::string::size_type end = g4ver.find(" ", start);
- g4ver = g4ver.substr(start, end - start);
- std::string g4major = g4ver.substr(1, 1);
- std::string g4minor = g4ver.substr(4, 1);
- std::string g4patch = "0";
- if (g4ver.find("patch") != string::npos) {
- g4patch = g4ver.substr(13, 1);
- }
- g4ver = g4major + g4minor + g4patch;
-#endif
- std::stringstream s;
- s << "geant4-" << "v" << g4ver[0] << "r" << g4ver[1] << "p" << g4ver[2];
- m_geant4VersionString = new std::string(s.str());
- }
- return *m_geant4VersionString;
-}
-
} // namespace slic
projects/slic/trunk/src
--- projects/slic/trunk/src/StdHepEventSource.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/StdHepEventSource.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,85 +1,66 @@
// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/StdHepEventSource.cc,v 1.10 2013-11-06 00:23:37 jeremy Exp $
#include "StdHepEventSource.hh"
-// slic
+// SLIC
#include "LcioManager.hh"
-#include "LcioMcpManager.hh"
+#include "MCParticleManager.hh"
+// LCIO
+#include "UTIL/LCTOOLS.h"
+
namespace slic {
StdHepEventSource::StdHepEventSource(const std::string& fname) :
- EventSourceWithInputFile("StdHepEventSource", fname), m_convertor(0) {
- // create loader with internal reader
- m_loader = new StdHepLoader();
-
- // create new convertor
- m_convertor = new StdHepToLcioConvertor();
+ EventSourceWithInputFile("StdHepEventSource", fname),
+ m_eventGenerator(0) {
}
void StdHepEventSource::generate(G4Event* anEvent) {
- LcioManager::instance()->getGenerator()->generatePrimaryVertexFromMcpCollection(
- LcioMcpManager::instance()->getInitialMcpCollection(), anEvent);
+ //std::cout << "StdHepEventSource::generate" << std::endl;
+
+ m_eventGenerator->generateEvent(anEvent);
}
// open the current file
void StdHepEventSource::open() {
- m_loader->openStdHepFile(m_filename);
+
+ /* 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();
+ delete m_eventGenerator;
m_fileIsOpen = false;
}
// read the next event
void StdHepEventSource::readNextEvent() {
- m_loader->readNextEvent();
- if (m_loader->isEndOfInput()) {
- m_eof = true;
+ 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();
-
- if (rdr) {
- rdr->printEventTable();
- } else {
- log() << LOG::error << "No current StdHep reader." << LOG::done;
- }
-
- log() << LOG::okay << LOG::endl;
+void StdHepEventSource::printCurrentEvent() {
+ if (MCParticleManager::instance()->getMCParticleCollection() != NULL)
+ UTIL::LCTOOLS::printMCParticles(MCParticleManager::instance()->getMCParticleCollection());
}
void StdHepEventSource::beginRun(const G4Run* aRun) {
+
// do superclass setup
EventSourceWithInputFile::beginRun(aRun);
-
- // setup convertor
- assert( m_loader->getStdHepReader());
- m_convertor->setStdHepReader(m_loader->getStdHepReader());
}
void StdHepEventSource::beginEvent(const G4Event* anEvent) {
- // read an event
- EventSourceWithInputFile::beginEvent(anEvent);
- // not EOF?
- if (!isEOF()) {
+ //std::cout << "StdHepEventSource::beginEvent" << std::endl;
- // set convertor's coll
- m_convertor->setParticleCollection(LcioMcpManager::instance()->getInitialMcpCollection());
-
- // run convertor
- m_convertor->fillParticleCollection();
- }
+ // read an event
+ EventSourceWithInputFile::beginEvent(anEvent);
}
-lStdHep* StdHepEventSource::getStdHepReader() {
- return m_loader->getStdHepReader();
}
-}
projects/slic/trunk/src
--- projects/slic/trunk/src/StdHepLoader.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/StdHepLoader.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,87 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/StdHepLoader.cc,v 1.23 2013-06-26 01:57:02 jeremy Exp $
-
-#include "StdHepLoader.hh"
-
-// lcdd
-#include "lcdd/util/StringUtil.hh"
-
-// geant4
-#include "G4PrimaryParticle.hh"
-#include "G4Event.hh"
-
-// std
-#include <assert.h>
-
-// lcio
-#include "UTIL/lStdHep.hh"
-
-using UTIL::lStdHep;
-
-namespace slic {
-
-StdHepLoader::~StdHepLoader() {
- ;
-}
-
-StdHepLoader::StdHepLoader() :
- Module("StdHepLoader"), m_filename(""), m_reader(0), m_isEndOfInput(false) {
- ;
-}
-
-void StdHepLoader::openStdHepFile(const std::string& filename) {
- openStdHepFile(filename.c_str());
-}
-
-void StdHepLoader::openStdHepFile(const char* filename) {
- // set filename
- m_filename = filename;
-
- // close if currently open
- closeStdHepFile();
-
- // new lStdHep interface
- m_reader = new lStdHep(filename);
-
- // if error, assume fatal
- if (m_reader->getError()) {
- G4Exception("StdHepLoader::openStdHepFile()", "", FatalException, "Error opening input file.");
- }
-
- // print open message
- log().okay(
- "StdHep file <" + std::string(filename) + "> opened with <"
- + StringUtil::toString((int) m_reader->numEvents()) + "> events.");
-}
-
-void StdHepLoader::closeStdHepFile() {
- if (0 != m_reader) {
- delete m_reader;
- m_reader = 0;
- }
-}
-
-void StdHepLoader::readNextEvent() {
- /* Read the next event and check the error code. */
- if (m_reader->readEvent() != LSH_SUCCESS) {
- /* The error code was not EOF. Assume this means a fatal error occurred. */
- if (m_reader->getError() != LSH_ENDOFFILE) {
- log().fatal("StdHep errorcode <" + StringUtil::toString((int) m_reader->getError()) + ">.");
- G4Exception("", "", FatalException, "lStdHep reported an unrecoverable error.");
- }
- /* EOF */
- else {
- m_isEndOfInput = true;
- }
- }
-}
-
-void StdHepLoader::dumpCurrentEvent() {
- if (m_reader) {
- m_reader->printEventHeader();
- m_reader->printEvent();
- m_reader->printEventTable();
- } else {
- log().error("Cannot run StdHepLoader::dumpCurrentEvent. StdHep reader is not initialized.");
- }
-}
-}
projects/slic/trunk/src
--- projects/slic/trunk/src/StdHepToLcioConvertor.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/StdHepToLcioConvertor.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,541 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/StdHepToLcioConvertor.cc,v 1.44 2013-11-06 00:23:37 jeremy Exp $
-#include "StdHepToLcioConvertor.hh"
-
-// slic
-#include "LcioMcpManager.hh"
-#include "HepPDTManager.hh"
-
-// lcio
-#include "IMPL/LCEventImpl.h"
-#include "UTIL/LCTOOLS.h"
-
-// geant4
-#include "G4ParticleDefinition.hh"
-#include "G4ParticleTable.hh"
-#include "G4PhysicalConstants.hh"
-
-// stl
-#include <cassert>
-#include <cmath>
-
-// slic
-#include "LcioMcpManager.hh"
-
-using EVENT::LCIO;
-using EVENT::LCCollection;
-using IMPL::LCEventImpl;
-using IMPL::MCParticleImpl;
-using UTIL::LCTOOLS;
-
-namespace slic {
-StdHepToLcioConvertor::StdHepToLcioConvertor(lStdHep* rdr) :
- Module("StdHepToLcioConverter"), m_reader(rdr) {
- assert(rdr);
-}
-
-StdHepToLcioConvertor::StdHepToLcioConvertor() :
- Module("StdHepToLcioConverter"), m_reader(0) {
-}
-
-StdHepToLcioConvertor::~StdHepToLcioConvertor() {
-}
-
-void StdHepToLcioConvertor::setParticleCollection(EVENT::LCCollection* mcpColl) {
- m_currentMcpColl = mcpColl;
-}
-
-EVENT::LCCollection* StdHepToLcioConvertor::getParticleCollection() {
- return m_currentMcpColl;
-}
-
-void StdHepToLcioConvertor::fillParticleCollection() {
- // Reset the MCParticle list.
- //m_currentMcpColl->clear();
- //G4cout << "m_currentMcpColl.size =" << m_currentMcpColl->getNumberOfElements() << G4endl;
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "********** Creating LCIO MCParticle Collection from StdHep **********" << LOG::done;
-#endif
-
- //assert( m_reader );
- //assert( m_currentMcpColl );
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "******** loop 1: create MCParticles ********" << LOG::done;
-#endif
-
- /* Loop #1 : Create MCParticles with no parentage or daughter links. */
- int ntracks = m_reader->nTracks();
- //G4cout << "StdHep.ntracks = " << ntracks << G4endl;
- for (int ihep = 0; ihep < ntracks; ihep++) {
-
-#ifdef SLIC_DEBUG
- printTrack( ihep );
-#endif
-
-#ifdef SLIC_DEBUG
- checkParentage( ihep );
-#endif
-
- // create MCP and add to coll
- createParticle(ihep);
-
-#ifdef SLIC_LOG
- log() << LOG::debug << LOG::endl << "--" << LOG::endl << LOG::done;
-#endif
- }
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "******** loop 2: parent setup ********" << LOG::done << LOG::done;
-#endif
-
- /* Loop #2 : Set the MCParticle parentage. */
- for (int ihep = 0; ihep < ntracks; ihep++) {
-
-#ifdef SLIC_LOG
- printIndex( ihep );
-// printMothers( ihep );
-#endif
-
- MCParticleImpl* mcp = dynamic_cast<MCParticleImpl*>(m_currentMcpColl->getElementAt(ihep));
-
- if (mcp) {
- setupParents(ihep, mcp);
- }
-
-#ifdef SLIC_LOG
- log() << LOG::debug << LOG::done << "--" << LOG::done << LOG::done;
-#endif
-
- }
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "******** loop 3: daughter setup ********" << LOG::done << LOG::done;
-#endif
-
- /* Loop #3 : Setup the daughters. */
- for (int ihep = 0; ihep < ntracks; ihep++) {
-
-#ifdef SLIC_LOG
- printIndex( ihep );
- printDaughters( ihep );
-#endif
-
- MCParticleImpl* mcp = dynamic_cast<MCParticleImpl*>(m_currentMcpColl->getElementAt(ihep));
-
- if (mcp) {
- setupDaughters(ihep, mcp);
- } else {
- log() << LOG::error << "mcp is null; ihep <" << ihep << "> does not appear to be a valid idx!"
- << LOG::done;
- }
-
-#ifdef SLIC_LOG
- log() << LOG::debug << LOG::done << "--" << LOG::done << LOG::done;
-#endif
-
- }
-
-#ifdef SLIC_LOG
- log() << LOG::debug << LOG::done;
-#endif
-}
-
-IMPL::MCParticleImpl* StdHepToLcioConvertor::createParticle(int ihep) {
- // new MCP
- MCParticleImpl* mcp = new MCParticleImpl();
-
- // PDG
- int pdgid = m_reader->pid(ihep);
- mcp->setPDG(pdgid);
-
- G4ParticleDefinition* pdef = G4ParticleTable::GetParticleTable()->FindParticle(pdgid);
-
-#ifdef SLIC_LOG
- if ( pdef != 0 ) {
- log() << LOG::debug << "found definition for particle <" << pdef->GetParticleName() << "> with pdgid <" << pdgid << ">" << LOG::done;
- }
- else {
- log() << LOG::debug << "no definition found for particle with PDG ID <" << pdgid << ">" << LOG::done;
- }
-#endif
-
- /*
- * Geant4 Particle definition exists for this MCParticle,
- * so set the charge from it.
- */
- if (pdef != 0) {
- mcp->setCharge(pdef->GetPDGCharge());
- }
- /*
- * No Geant4 particle definition, so flag charge as invalid.
- */
- else {
- mcp->setCharge(LcioMcpManager::m_NAN);
- }
-
- // momentum vec
- float p[3] = { m_reader->Px(ihep), m_reader->Py(ihep), m_reader->Pz(ihep) };
- mcp->setMomentum(p);
-
- // mass
- mcp->setMass(m_reader->M(ihep));
-
- // vertex
- double vtx[3] = { m_reader->X(ihep), m_reader->Y(ihep), m_reader->Z(ihep) };
- mcp->setVertex(vtx);
-
- // generator status
- mcp->setGeneratorStatus(m_reader->status(ihep));
-
- // sim status = none
- mcp->setSimulatorStatus(0);
-
- // creation time
- mcp->setTime(m_reader->T(ihep) / c_light);
-
-#if LCIO_VERSION_GE(1, 60)
-
- // Spin and color flow code only applicable to StdHep4 files.
- if(m_reader->isStdHepEv4())
- {
- // Set particle spin.
- float spin[3];
- spin[0] = m_reader->spinX(ihep);
- spin[1] = m_reader->spinY(ihep);
- spin[2] = m_reader->spinZ(ihep);
- mcp->setSpin(spin);
-
- // Set particle color flow.
- int colorFlow[2] = {m_reader->colorflow(ihep, 0), m_reader->colorflow(ihep, 1)};
- mcp->setColorFlow(colorFlow);
- }
-#endif
-
- // Add particle to collection.
- m_currentMcpColl->addElement(mcp);
-
- return mcp;
-}
-
-// util functions
-void StdHepToLcioConvertor::setupParents(int ihep, IMPL::MCParticleImpl* mcp) {
-//#ifdef SLIC_LOG
-// log() << LOG::debug << "setupParents" << LOG::done;
-//#endif
-
- // get parent indices
- int mom1_idx = m_reader->mother1(ihep) - 1;
- int mom2_idx = m_reader->mother2(ihep) - 1;
-
- std::vector<int> parVec;
- makeIndexVec(mom1_idx, mom2_idx, parVec);
-
- if (parVec.size() > 0) {
- for (std::vector<int>::iterator it = parVec.begin(); it != parVec.end(); it++) {
- if (*it != ihep) {
- addParent(*it, mcp);
- }
-//#ifdef SLIC_LOG
-// else
-// {
-// log() << LOG::debug << "Ignoring parent pointing to self for ihep <" << ihep << ">" << LOG::done;
-// }
-//#endif
- }
- }
-}
-
-void StdHepToLcioConvertor::setupDaughters(int ihep, IMPL::MCParticleImpl* mcp) {
-//#ifdef SLIC_LOG
-// log() << LOG::debug << "setupDaughters for ihep <" << ihep << "> and MCParticle <" << mcp << ">" << LOG::done;
-//#endif
-#ifdef SLIC_LOG
- log() << LOG::debug << "MCParticle idx: " << LcioMcpManager::instance()->getParticleIndex(LcioMcpManager::instance()->getInitialMcpCollection(), mcp) << LOG::done;
-#endif
-
- // get dau indices
- int dau1_idx = m_reader->daughter1(ihep) % 10000 - 1;
- int dau2_idx = m_reader->daughter2(ihep) % 10000 - 1;
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "dau1 <" << dau1_idx << ">" << LOG::done;
- log() << LOG::debug << "dau2 <" << dau2_idx << ">" << LOG::done;
-#endif
-
- // Check for valid daughter indices. Some generators such as Whizard put bad values in these fields.
- if (!isValidIndex(dau1_idx) || !isValidIndex(dau2_idx)) {
-#ifdef SLIC_LOG
- log() << LOG::debug << "ihep <" << ihep << "> has invalid daughter index range <" << dau1_idx << ", " << dau2_idx << ">." << LOG::done;
-#endif
- return;
- }
-
- std::vector<int> dauVec;
- dauVec = makeIndexVec(dau1_idx, dau2_idx, dauVec);
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "printing dau vec ..." << LOG::done;
- for (std::vector<int>::const_iterator it = dauVec.begin();
- it != dauVec.end();
- it++)
- {
- log() << LOG::debug << (*it) << " ";
- }
- log () << LOG::done;
-
-#endif
-
- if (dauVec.size() > 0) {
- for (std::vector<int>::iterator it = dauVec.begin(); it != dauVec.end(); it++) {
- if (*it != ihep) {
- addDaughter(*it, mcp);
- }
-#ifdef SLIC_LOG
- else
- {
- log() << LOG::error << "Ignoring daughter pointing to self!" << LOG::done;
- }
-#endif
- }
- }
-#ifdef SLIC_LOG
- else {
- log() << LOG::debug << "particle has no daughters" << LOG::done;
- }
-#endif
-}
-
-bool StdHepToLcioConvertor::hasParent(IMPL::MCParticleImpl* dauMcp, IMPL::MCParticleImpl* parMcp) {
- bool isPar = false;
- if (dauMcp && parMcp) {
- int numParents = dauMcp->getParents().size();
-
- MCParticleImpl* pMomMcp;
-
- // find if already lists as parent
- for (int i = 0; i < numParents; i++) {
-
- pMomMcp = static_cast<MCParticleImpl*>(dauMcp->getParents()[i]);
-
- if (pMomMcp == parMcp) {
- isPar = true;
- break;
- }
- }
- }
-
- return isPar;
-}
-
-void StdHepToLcioConvertor::addDaughter(int dauIdx, IMPL::MCParticleImpl* parMcp) {
- if (dauIdx >= 0) {
-
- MCParticleImpl* dauMcp = static_cast<MCParticleImpl*>(m_currentMcpColl->getElementAt(dauIdx));
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "found dau particle <" << dauMcp << ">" << LOG::done;
-#endif
-
-// assert( dauMcp );
-
- if (dauMcp) {
- if (!hasParent(dauMcp, parMcp)) {
- dauMcp->addParent(parMcp);
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "added daughter <" << dauIdx << ">" << LOG::done;
-#endif
- }
-#ifdef SLIC_LOG
- else {
- log() << LOG::debug << "daughter <" << dauIdx << "> already has this parent" << LOG::done;
- }
-#endif
- }
-#ifdef SLIC_LOG
- else {
- log() << LOG::error << "WARNING: dauMcp or parMcp is null!" << LOG::done;
- }
-#endif
-
- }
-}
-
-void StdHepToLcioConvertor::addParent(int parIdx, IMPL::MCParticleImpl* mcp) {
-#ifdef SLIC_LOG
- log() << LOG::debug << "addMcpParent" << LOG::done;
-#endif
-
- if (parIdx >= 0) {
-
- /* If index is > size of collection, the particle's parent doesn't exist yet. We need to die! */
- if (parIdx > m_currentMcpColl->getNumberOfElements()) {
-#ifdef SLIC_LOG
- log() << LOG::fatal << "StdHepToLcioConverter::addMcpParent - Parent index is out of range <" << parIdx << ">" << LOG::done;
-#endif
- G4Exception("", "", FatalException, "StdHep index (ihep) is out of range.");
- }
-
- MCParticleImpl* parMcp = static_cast<MCParticleImpl*>(m_currentMcpColl->getElementAt(parIdx));
-
- if (parMcp) {
- if (!hasParent(mcp, parMcp)) {
- mcp->addParent(parMcp);
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "added parent <" << parIdx << ">" << LOG::done;
-#endif
- }
-#ifdef SLIC_LOG
- else {
- log() << LOG::debug << "daughter already has parent <" << parIdx << ">" << LOG::done;
- }
-#endif
- } else {
- G4Exception("", "", FatalException, "Failed to get parent particle from MCParticle collection.");
- }
- }
-#ifdef SLIC_LOG
- else {
- log() << LOG::debug << "ignoring parIdx = -1 " << LOG::done;
- }
-#endif
-}
-
-std::vector<int> StdHepToLcioConvertor::makeIndexVec(int idx1, int idx2, std::vector<int>& vec) {
-#ifdef SLIC_LOG
- log() << LOG::debug << "idx1 <" << idx1 << ">" << LOG::done;
- log() << LOG::debug << "idx2 <" << idx2 << ">" << LOG::done;
-#endif
-
- if (idx1 >= 0 && idx2 >= 0) {
-
- if (idx1 < idx2) {
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "range: idx1 to idx2" << LOG::done;
-#endif
-
- for (int i = idx1; i < (idx2 + 1); i++) {
- vec.push_back(i);
- }
- } else if (idx1 > idx2) {
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "discrete: idx1 and idx2" << LOG::done;
-#endif
-
- vec.push_back(idx1);
- vec.push_back(idx2);
- }
- // indices are equal
- else {
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "single: idx1 == idx2" << LOG::done;
-#endif
-
- vec.push_back(idx1);
- }
- } else if (idx1 >= 0) {
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "single: idx1 only" << LOG::done;
-#endif
-
- vec.push_back(idx1);
- } else if (idx2 >= 0) {
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "single: idx2 only" << LOG::done;
-#endif
-
- vec.push_back(idx2);
- }
-
- return vec;
-}
-
-void StdHepToLcioConvertor::printIndex(int ihep) {
- log() << LOG::debug << "ihep <" << ihep << ">" << LOG::done;
-}
-
-void StdHepToLcioConvertor::printMothers(int ihep) {
- log() << LOG::debug << "mom1 <" << m_reader->mother1(ihep) - 1 << ">" << LOG::done;
- log() << LOG::debug << "mom2 <" << m_reader->mother2(ihep) - 1 << ">" << LOG::done;
-}
-
-void StdHepToLcioConvertor::printDaughters(int ihep) {
- log() << LOG::debug << "dau1 <" << m_reader->daughter1(ihep) - 1 << ">" << LOG::done;
- log() << LOG::debug << "dau2 <" << m_reader->daughter2(ihep) - 1 << ">" << LOG::done;
-}
-
-void StdHepToLcioConvertor::printTrack(int ihep) {
- printIndex(ihep);
-
- log() << LOG::debug << "pid <" << m_reader->pid(ihep) << ">" << LOG::done;
- log() << LOG::debug << "M = " << m_reader->M(ihep) << LOG::done;
- log() << LOG::debug << "T = " << m_reader->T(ihep) << LOG::done;
- log() << LOG::debug << "status <" << m_reader->status(ihep) << ">" << LOG::done;
-
- printMothers(ihep);
- printDaughters(ihep);
-
- log() << LOG::debug << "P = (" << m_reader->Px(ihep) << ", " << m_reader->Py(ihep) << ", "
- << m_reader->Pz(ihep) << ")" << LOG::done;
-
- log() << LOG::debug << "Vtx = (" << m_reader->X(ihep) << ", " << m_reader->Y(ihep) << ", "
- << m_reader->Z(ihep) << ")" << LOG::done;
-
- log() << LOG::debug << LOG::done;
-}
-
-void StdHepToLcioConvertor::checkParentage(int ihep) {
- int mom1 = m_reader->mother1(ihep) - 1;
- int mom2 = m_reader->mother2(ihep) - 1;
- int dau1 = m_reader->daughter1(ihep) - 1;
- int dau2 = m_reader->daughter2(ihep) - 1;
-
- // dau1 < ihep
- if (dau1 >= 0 && dau1 < ihep) {
- log() << LOG::debug << "WARNING: dau1 < ihep" << LOG::done;
- }
-
- // dau2 < ihep
- if (dau2 >= 0 && dau2 < ihep) {
- log() << LOG::debug << "WARNING: dau2 < ihep" << LOG::done;
- }
-
- // mom1 > ihep
- if (mom1 >= 0 && mom1 > ihep) {
- log() << LOG::debug << "WARNING: mom1 > ihep" << LOG::done;
- }
-
- // mom2 > ihep
- if (mom2 >= 0 && mom2 > ihep) {
- log() << LOG::debug << "WARNING: mom2 > ihep" << LOG::done;
- }
-
- // first particle in list has parents
- if (ihep == 0 && (mom1 >= 0 || mom2 >= 0)) {
- log() << LOG::debug << "WARNING: ihep == 0 has parents" << LOG::done;
- }
-
- // particle past first 2 with no mother
- if (ihep > 1 && (mom1 < 0 && mom2 < 0)) {
- log() << LOG::debug << "WARNING: ihep > 1 with no mother" << LOG::done;
- }
-
- // mother points to self
- if (ihep == mom1 || ihep == mom2) {
- log() << LOG::debug << "WARNING: mom1 or mom2 == ihep" << LOG::done;
- }
-
- // daughter points to self
- if (ihep == dau1 || ihep == dau2) {
- log() << LOG::debug << "WARNING: dau1 or dau2 == ihep" << LOG::done;
- }
-}
-
-} // namespace slic
projects/slic/trunk/src
--- projects/slic/trunk/src/SteppingAction.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/SteppingAction.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -2,13 +2,14 @@
// SLIC
#include "SteppingAction.hh"
-#include "TrajectoryManager.hh"
+#include "UserTrackInformation.hh"
+#include "TrackUtil.hh"
// LCDD
#include "lcdd/geant4/UserRegionInformation.hh"
-#include "lcdd/hits/TrackInformation.hh"
// Geant4
+#include "G4SteppingManager.hh"
#include "G4UserSteppingAction.hh"
#include "G4Track.hh"
@@ -20,16 +21,38 @@
SteppingAction::~SteppingAction() {
}
-void SteppingAction::UserSteppingAction(const G4Step* aStep) {
- /* TrajectoryManager's stepping action */
- TrajectoryManager::instance()->stepping(aStep);
+void SteppingAction::UserSteppingAction(const G4Step* step) {
+ /* Check for back scattering. */
+ if (isBackScattering(step)) {
+
+ // Get track information.
+ UserTrackInformation* trackInfo =
+ (UserTrackInformation*)fpSteppingManager->GetTrack()->GetUserInformation();
+
+ // Need to start making track information if doesn't exist because these tracks or daughters could leave tracker hits.
+ if (trackInfo == NULL)
+ // Add new track information but don't save. The tracker hit processor can turn saving on.
+ trackInfo = TrackUtil::setupUserTrackInformation(fpSteppingManager->GetTrack(), false);
+
+ TrackSummary* trackSummary = trackInfo->getTrackSummary();
+
+ //G4cout << "setting backscattering flag" << G4endl;
+ trackInfo->getTrackSummary()->update(fpSteppingManager->GetTrack());
+ trackInfo->getTrackSummary()->setBackScattering();
+
+ // Suspend backscattering tracks to postpone them to the end of shower development.
+ fpSteppingManager->GetTrack()->SetTrackStatus(fSuspend);
+
+ //G4cout << "suspended backscattering track " << fpSteppingManager->GetTrack()->GetTrackID() << G4endl;
+ }
+
/* Check if the track should be killed. */
- checkKillTrack(aStep);
+ checkKillTrack(step);
}
void SteppingAction::checkKillTrack(const G4Step* aStep) {
- /* Check if the preStepPoint region has tracking killing enabled. */
+ /* Check if the region of the pre-point has tracking killing enabled. */
UserRegionInformation* regionInfo = UserRegionInformation::getRegionInformation(aStep->GetPreStepPoint());
if (regionInfo->getKillTracks()) {
/* Kill the track if region track kill flag is set to true. */
@@ -38,4 +61,25 @@
}
}
+bool SteppingAction::isBackScattering(const G4Step* step) {
+
+ /* Get pre-point info. */
+ const G4StepPoint* prePoint = step->GetPreStepPoint();
+ const G4StepPoint* postPoint = step->GetPostStepPoint();
+ const UserRegionInformation* prePointRegionInfo =
+ (UserRegionInformation*)prePoint->GetTouchableHandle()->GetVolume()->GetLogicalVolume()->GetRegion()->GetUserInformation();
+ G4bool prePointInsideTrackerRegion = prePointRegionInfo->getStoreSecondaries();
+
+ /* Get post-point info. */
+ G4bool postPointInsideTrackerRegion = false;
+ if (postPoint->GetTouchableHandle()->GetVolume() != NULL) {
+ const UserRegionInformation* postPointRegionInfo =
+ (UserRegionInformation*)postPoint->GetTouchableHandle()->GetVolume()->GetLogicalVolume()->GetRegion()->GetUserInformation();
+ postPointInsideTrackerRegion = postPointRegionInfo->getStoreSecondaries();
+ }
+
+ /* Track is backscattering if it steps from a non-tracking region back into a tracking region. */
+ return ((!prePointInsideTrackerRegion) && postPointInsideTrackerRegion);
+}
+
} // namespace slic
projects/slic/trunk/src
--- projects/slic/trunk/src/TrackingAction.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/TrackingAction.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,28 +1,91 @@
// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/TrackingAction.cc,v 1.14 2012-11-27 19:32:19 jeremy Exp $
#include "TrackingAction.hh"
-// slic
-#include "Trajectory.hh"
-#include "TrajectoryManager.hh"
+// LCDD
+#include "lcdd/detectors/CurrentTrackState.hh"
+// SLIC
+#include "TrackUtil.hh"
+#include "UserTrackInformation.hh"
+
+// Geant4
+#include "G4TransportationManager.hh"
+
namespace slic {
-TrackingAction::TrackingAction() {
- // set TrajectoryManager links
- m_trjMgr = TrajectoryManager::instance();
-}
+void TrackingAction::PreUserTrackingAction(const G4Track* track) {
-TrackingAction::~TrackingAction() {
+ G4int trackID = track->GetTrackID();
+ G4int parentID = track->GetParentID();
+
+ G4bool isTruePrimary = (track->GetParentID() == 0);
+ UserRegionInformation* regionInfo = TrackUtil::getRegionInformation(track);
+ G4bool insideTrackerRegion = regionInfo->getStoreSecondaries();
+ UserTrackInformation* trackInfo = TrackUtil::getUserTrackInformation(track);
+ bool hasTrackInfo = (trackInfo != NULL);
+ bool aboveEnergyThreshold = (track->GetKineticEnergy() > regionInfo->getThreshold());
+
+ // This happens when backscattering track are un-suspended.
+ if (hasTrackInfo) {
+ // Track is backscattering?
+ if (trackInfo->getTrackSummary()->getBackScattering()) {
+ // Update current track ID for calorimeter SDs and return.
+ CurrentTrackState::setCurrentTrackID(trackID);
+ return;
+ }
+ return;
+ }
+
+ // Primary OR in tracker region and above energy threshold?
+ if (isTruePrimary || (insideTrackerRegion && aboveEnergyThreshold)) {
+
+ // Update track ID for calorimeter SDs.
+ CurrentTrackState::setCurrentTrackID(trackID);
+
+ // Turn on trajectory storage.
+ fpTrackingManager->SetStoreTrajectory(true);
+
+ // Has track information?
+ if (!hasTrackInfo) {
+ // Create the track information and save it.
+ TrackUtil::setupUserTrackInformation(track, true);
+ }
+
+ } else {
+ // Inside tracker region and below energy threshold?
+ if (insideTrackerRegion && !aboveEnergyThreshold) {
+ // Has track information?
+ if (!hasTrackInfo) {
+ // Create track information but do not save it.
+ TrackUtil::setupUserTrackInformation(track, false);
+ }
+ } else {
+ // Track is outside tracking region so turn off trajectory storage.
+ // These are generally tracks from calorimeter showers.
+ fpTrackingManager->SetStoreTrajectory(false);
+ }
+ }
}
-void TrackingAction::PreUserTrackingAction(const G4Track* aTrack) {
- // TrajectoryManager's action; also requires setup of trkMgr
- m_trjMgr->setTrackingManager(fpTrackingManager); // needs to be called every time?
- m_trjMgr->preTracking(aTrack);
+void TrackingAction::PostUserTrackingAction(const G4Track* track) {
+
+ // Get the track information and summary.
+ UserTrackInformation* trackInfo = TrackUtil::getUserTrackInformation(track);
+ TrackSummary* trackSummary = NULL;
+ if (trackInfo != NULL)
+ trackSummary = trackInfo->getTrackSummary();
+
+ // Track has secondaries?
+ if (fpTrackingManager->GimmeSecondaries()->size() > 0 && trackInfo == NULL) {
+ // Tracks with secondaries always get track information but they are not saved by default.
+ trackSummary = TrackUtil::setupUserTrackInformation(track, false)->getTrackSummary();
+ }
+
+ // Track summary exists?
+ if (trackSummary != NULL) {
+ // Update the track summary.
+ trackSummary->update(track);
+ }
}
-void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) {
- // TrajectoryManager's action
- m_trjMgr->postTracking(aTrack);
-}
-}
+} // namespace slic
projects/slic/trunk/src
--- projects/slic/trunk/src/Trajectory.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/Trajectory.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,127 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/Trajectory.cc,v 1.17 2013-06-26 01:57:02 jeremy Exp $
-#include "Trajectory.hh"
-
-// lcdd
-#include "lcdd/hits/TrackInformation.hh"
-
-// geant4
-#include "G4TrajectoryPoint.hh"
-
-namespace slic {
-
-G4Allocator<Trajectory> TrajectoryAllocator;
-
-Trajectory::Trajectory() :
- G4VTrajectory(), m_particleDefinition(0), m_particleName(""), m_PDGCharge(0), m_PDGEncoding(0), m_trackID(
- 0), m_parentID(0), m_posRecord(0), m_globalTime(0), m_endpointEnergy(0), m_createdInSimulation(
- false), m_hasEndpointDaughters(false), m_backscatter(false), m_vertexIsNotEndpointOfParent(
- false) {
-}
-
-Trajectory::Trajectory(const G4Track *aTrack) :
- G4VTrajectory() {
- m_particleDefinition = aTrack->GetDefinition();
- m_particleName = m_particleDefinition->GetParticleName();
- m_PDGCharge = m_particleDefinition->GetPDGCharge();
- m_PDGEncoding = m_particleDefinition->GetPDGEncoding();
- m_trackID = aTrack->GetTrackID();
- m_parentID = aTrack->GetParentID();
- m_posRecord = new TrajectoryPointContainer();
- m_posRecord->push_back(new G4TrajectoryPoint(aTrack->GetPosition()));
- m_initialMomentum = aTrack->GetMomentum();
-
- // default flags to false
- m_createdInSimulation = false;
- m_hasEndpointDaughters = false;
- m_backscatter = false;
- m_vertexIsNotEndpointOfParent = false;
-
- m_globalTime = aTrack->GetGlobalTime();
-}
-
-void Trajectory::setupManuallyFromTrack(const G4Track* aTrack) {
- TrackInformation* trkInfo = static_cast<TrackInformation*>(aTrack->GetUserInformation());
-
- assert( trkInfo);
-
- m_parentID = trkInfo->getOriginalTrackID();
-
- // position already added in trajectory ctor
- m_posRecord->insert(m_posRecord->begin(), new G4TrajectoryPoint(aTrack->GetVertexPosition()));
-
- G4double vKE = aTrack->GetVertexKineticEnergy();
-
- G4double mass = m_particleDefinition->GetPDGMass();
- G4double magmom = sqrt(vKE * vKE + 2. * vKE * mass);
- m_initialMomentum = magmom * (aTrack->GetVertexMomentumDirection());
-
- m_globalTime = aTrack->GetGlobalTime();
-}
-
-Trajectory::Trajectory(Trajectory &right) :
- G4VTrajectory() {
- m_particleName = right.m_particleName;
- m_particleDefinition = right.m_particleDefinition;
- m_PDGCharge = right.m_PDGCharge;
- m_PDGEncoding = right.m_PDGEncoding;
- m_trackID = right.m_trackID;
- m_parentID = right.m_parentID;
- m_globalTime = right.m_globalTime;
-
- m_posRecord = new TrajectoryPointContainer();
- for (unsigned int i = 0; i < right.m_posRecord->size(); i++) {
- G4TrajectoryPoint *rightPoint = (G4TrajectoryPoint*) ((*(right.m_posRecord))[i]);
- m_posRecord->push_back(new G4TrajectoryPoint(*rightPoint));
- }
-}
-
-Trajectory::~Trajectory() {
- for (unsigned int i = 0; i < m_posRecord->size(); i++) {
- delete (*m_posRecord)[i];
- }
-
- m_posRecord->clear();
- delete m_posRecord;
-}
-
-void Trajectory::AppendStep(const G4Step *aStep) {
- m_posRecord->push_back(new G4TrajectoryPoint(aStep->GetPostStepPoint()->GetPosition()));
-}
-
-void Trajectory::MergeTrajectory(G4VTrajectory *secondTrajectory) {
- if (!secondTrajectory)
- return;
-
- Trajectory *seco = (Trajectory*) secondTrajectory;
- G4int pent = seco->GetPointEntries();
- for (int i = 1; i < pent; i++) {
- m_posRecord->push_back((*(seco->m_posRecord))[i]);
- }
-
- delete (*seco->m_posRecord)[0];
- seco->m_posRecord->clear();
-}
-
-const std::string& Trajectory::getFinalStatusString(EFinalStatus fs) {
- static const std::string str_unset = "Unset";
- static const std::string str_invalid = "Invalid";
- static const std::string str_decayed_trk = "DecayedInTracker";
- static const std::string str_left = "LeftDetector";
- static const std::string str_stopped = "Stopped";
- static const std::string str_decayed_cal = "DecayedInCalorimeter";
-
- if (fs == eDecayedInTracker) {
- return str_decayed_trk;
- } else if (fs == eLeftDetector) {
- return str_left;
- } else if (fs == eStopped) {
- return str_stopped;
- } else if (fs == eDecayedInCalorimeter) {
- return str_decayed_cal;
- } else if (fs == eUnset) {
- return str_unset;
- }
-
- return str_invalid;
-}
-}
projects/slic/trunk/src
--- projects/slic/trunk/src/TrajectoryManager.cc 2014-07-15 23:45:13 UTC (rev 3186)
+++ projects/slic/trunk/src/TrajectoryManager.cc 2014-07-15 23:59:54 UTC (rev 3187)
@@ -1,584 +0,0 @@
-// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/TrajectoryManager.cc,v 1.54 2013-11-06 00:23:37 jeremy Exp $
-
-// slic
-#include "TrajectoryManager.hh"
-#include "Trajectory.hh"
-
-// geant4
-#include "G4TrackVector.hh"
-#include "G4Event.hh"
-
-namespace slic {
-
-TrajectoryManager::~TrajectoryManager() {
-}
-
-TrajectoryManager::TrajectoryManager() :
- Module("TrajectoryManager"), m_trackingManager(0), m_currentTrajectory(0), m_currentTrack(0), m_currentTrackID(
- -1) {
-}
-
-void TrajectoryManager::preTracking(const G4Track* aTrack) {
- m_currentTrack = aTrack;
- m_currentTrackID = aTrack->GetTrackID();
-
-#ifdef SLIC_LOG
- log() << LOG::verbose
- << "TrajectoryManager::preTracking()" << LOG::endl
- << "trackID <" << aTrack->GetTrackID() << ">" << LOG::endl
- << "particle <" << aTrack->GetDefinition()->GetParticleName() << ">" << LOG::endl
- << "volume <" << aTrack->GetVolume()->GetName() << ">"
- << LOG::done;
- printStatusFlags(aTrack);
-#endif
-
- // default trajectory storage to OFF
- m_trackingManager->SetStoreTrajectory(false);
-
- // handle primary or secondary track
- if (isPrimaryTrack(aTrack)) {
- handlePrimaryTrack(aTrack);
- } else {
- handleSecondaryTrack(aTrack);
- }
-}
-
-void TrajectoryManager::printStatusFlags(const G4Track* trk) {
- log() << LOG::debug << "status flags for trackID <" << trk->GetTrackID() << ">" << LOG::done;
- printStatusFlags(TrackInformation::getTrackInformation(trk));
-}
-
-void TrajectoryManager::printStatusFlags(TrackInformation* trkInfo) {
- if (trkInfo) {
- log() << LOG::verbose << "tracking status <"
- << TrackInformation::getTrackingStatusString(trkInfo->getTrackingStatus()) << ">" << LOG::endl
- << "orig tracking status <"
- << TrackInformation::getTrackingStatusString(trkInfo->getOriginalTrackingStatus()) << ">"
- << LOG::endl << "backscatter flag <" << trkInfo->getBackscatter() << ">" << LOG::endl
- << "belowThreshold flag <" << trkInfo->getBelowThreshold() << ">" << LOG::endl
- << "tracker hit flag <" << trkInfo->hasTrackerHit() << ">" << LOG::endl << "vtxnotendp flag <"
- << trkInfo->getVertexIsNotEndpointOfParent() << ">" << LOG::done;
- }
-#ifdef SLIC_LOG
- else {
- log() << LOG::debug << "Track info NULL in printStatusFlags()." << LOG::done;
- }
-#endif
-}
-
-void TrajectoryManager::printTrajectory(const Trajectory* trj) {
- log() << LOG::debug << "trajectory info" << LOG::endl << "trackID <" << trj->GetTrackID() << ">"
- << LOG::endl << "final status: " << trj->getFinalStatusString() << LOG::endl << "endp dau: "
- << trj->getHasEndpointDaughters() << LOG::endl << "backscatter: " << trj->getBackscatter()
- << LOG::endl << "vtxnotendp: " << trj->getVertexIsNotEndpointOfParent() << LOG::endl
- << "created in sim: " << trj->getCreatedInSimulation() << LOG::endl << "vtx: "
- << trj->GetPoint(0)->GetPosition() << LOG::endl << "endp: "
- << trj->GetPoint(trj->GetPointEntries() - 1)->GetPosition() << LOG::endl << LOG::done;
-}
-
-void TrajectoryManager::postTracking(const G4Track* aTrack) {
-#ifdef SLIC_LOG
- log() << LOG::verbose << "TrajectoryManager::postTracking()" << LOG::endl;
-#endif
-
- // check if needs manual setup
- if (needsManualTrajectorySetup(aTrack)) {
-
-#ifdef SLIC_LOG
- log() << LOG::verbose << "needs manual trajectory setup" << LOG::done;
-#endif
-
- // setup trajectory in default way
- Trajectory* trj = createTrajectoryFromTrack(aTrack);
-
- // do additional, manual setup using this track
- trj->setupManuallyFromTrack(aTrack);
- }
-
- // trk info to secondaries
- copySecondaryTrackInformationFromParent(aTrack);
-
- // set final status codes in current trajectory
- setTrajectoryFinalStatus(aTrack);
-
- // print if store traj on/off
-#ifdef SLIC_LOG
- log() << LOG::verbose << "store trajectory <" << m_trackingManager->GetStoreTrajectory() << ">" << LOG::done;
-#endif
-
- // print status flags
-#ifdef SLIC_LOG
- printStatusFlags(aTrack);
-#endif
-
- Trajectory* trj = static_cast<Trajectory*>(m_trackingManager->GimmeTrajectory());
- if (0 != trj) {
- trj->setEndpointEnergy(aTrack->GetTotalEnergy());
- }
-}
-
-void TrajectoryManager::handlePrimaryTrack(const G4Track* aTrack) {
- // set track info
- setPrimaryTrackInformation(aTrack);
-
- // create the trajectory
- createTrajectoryFromTrack(aTrack);
-}
-
-Trajectory* TrajectoryManager::createTrajectoryFromTrack(const G4Track* aTrack) {
-#ifdef SLIC_LOG
- log() << LOG::verbose << "createTrajectoryFromTrack()" << LOG::done;
-#endif
- Trajectory* trj = new Trajectory(aTrack);
- m_trackingManager->SetStoreTrajectory(true);
- m_trackingManager->SetTrajectory(trj);
- addTrackIDToTrajectoryLink(trj);
- m_currentTrajectory = trj;
- return trj;
-}
-
-void TrajectoryManager::setPrimaryTrackInformation(const G4Track* aTrack) {
- TrackInformation* trkInfo = new TrackInformation(aTrack);
-
- // set tracking flag using reg info
- setTrackingFlagFromRegionInformation(trkInfo, UserRegionInformation::getRegionInformation(aTrack));
-
- G4Track* nconstTrack = const_cast<G4Track*>(aTrack);
- nconstTrack->SetUserInformation(trkInfo);
-}
-
-void TrajectoryManager::setTrackingFlagFromRegionInformation(TrackInformation* trk_info,
- UserRegionInformation* reg_info) {
- if (reg_info->getStoreSecondaries() == true) {
- trk_info->setTrackingStatus(TrackInformation::eInTrackingRegion);
- } else {
- trk_info->setTrackingStatus(TrackInformation::eInNontrackingRegion);
- }
-}
-
-/*
- Called in PreTracking to set status flags and
- decide whether or not to create a trj for
- this secondary particle.
- */
-void TrajectoryManager::handleSecondaryTrack(const G4Track* aTrack) {
- // get track information
- TrackInformation* trkInfo = TrackInformation::getTrackInformation(aTrack);
-
- // get region information
- UserRegionInformation* regInfo = UserRegionInformation::getRegionInformation(aTrack);
-
- // must have track and region info to proceed
- if (regInfo && trkInfo) {
-
- // check if in tracking region
- if (regInfo->getStoreSecondaries()) {
-
- // check if backscatter
- if (isBackscatter(aTrack)) {
- trkInfo->setBackscatter(true);
- }
-
- /*
- Check if below threshold and continuous process creation,
- i.e. not an endp dau.
- */
- else if (isBelowThreshold(aTrack) && trkInfo->getVertexIsNotEndpointOfParent()) {
- trkInfo->setBelowThreshold(true);
- }
-
- /*
- Process of elimination: above threshold and not backscatter,
- so requires a new trajectory and trkInfo.
- */
- else {
-
- // setup secondary track info
- setupSecondaryTrackInformation(aTrack);
-
- // create the secondary trajectory
- createSecondaryTrajectory(aTrack);
- }
-
- /*
- // DEBUG
- if ( trkInfo != TrackInformation::getTrackInformation( aTrack ) ) {
- log() << "WARNING: track info was reset" << LOG::endl;
-
- log() << LOG::endl;
- log() << "old flags" << LOG::endl;
- printStatusFlags( trkInfo );
- log() << LOG::endl;
-
- log() << "new flags" << LOG::endl;
- printStatusFlags( TrackInformation::getTrackInformation( aTrack ) );
- log() << LOG::endl;
- }
- */
-
- /*
- Set flag for in tracking region, re-getting track info,
- because may have been recreated above.
- */
- TrackInformation::getTrackInformation(aTrack)->setTrackingStatus(
- TrackInformation::eInTrackingRegion);
- }
- // in calorimeter, i.e. non-tracking
- else {
-
- // set non tracking status
- trkInfo->setTrackingStatus(TrackInformation::eInNontrackingRegion);
- trkInfo->setOriginalTrackingStatus(TrackInformation::eInNontrackingRegion);
- }
- }
- // fatal error: regInfo or trkInfo were null
- else {
- log() << LOG::fatal << "FATAL ERROR: region info or track info were null in handleSecondaryTrack()."
- << LOG::endl << "region info <" << regInfo << ">" << LOG::endl << "track info <" << trkInfo
- << ">" << LOG::done;
- G4Exception("", "", FatalException, "Track info or region info was null.");
- }
-}
-
-void TrajectoryManager::setupSecondaryTrackInformation(const G4Track* aTrack) {
- // old track info
- TrackInformation* oldTrkInfo = TrackInformation::getTrackInformation(aTrack);
-
- // create new track info
- TrackInformation* newTrkInfo = new TrackInformation(aTrack);
-
- // copy vtxnotendp flag from old trkInfo
- newTrkInfo->setVertexIsNotEndpointOfParent(oldTrkInfo->getVertexIsNotEndpointOfParent());
-
- /*
- Delete old track info to avoid mem leak.
- */
- delete oldTrkInfo;
- oldTrkInfo = 0;
-
- /*
- Set original tracking status; it must be tracking, because otherwise
- would not be calling this function. (By definition, secondary tracks
- do not get trajectories in non-tracking regions.)
- */
- newTrkInfo->setOriginalTrackingStatus(TrackInformation::eInTrackingRegion);
-
- /*
- Setup new track info in this track.
- */
- G4Track* nconstTrk = const_cast<G4Track*>(aTrack);
- nconstTrk->SetUserInformation(newTrkInfo);
-}
-
-Trajectory* TrajectoryManager::createSecondaryTrajectory(const G4Track* aTrack) {
- // create the trajectory
- Trajectory* trj = createTrajectoryFromTrack(aTrack);
-
- // set CreatedInSimulation flag
- trj->setCreatedInSimulation(true);
-
- return trj;
-}
-
-// called in PostTracking to pass trk info to the track's secondaries
-void TrajectoryManager::copySecondaryTrackInformationFromParent(const G4Track* aTrack) {
- // get track info for parent
- TrackInformation* parTrkInfo = TrackInformation::getTrackInformation(aTrack);
-
- // get current trajectory from tracking manager
- Trajectory* parTrj = getCurrentTrajectory();
-
- // get secondaries container
- G4TrackVector* secondaries = m_trackingManager->GimmeSecondaries();
-
- // loop over secondaries to 1) create trk info and 2) set vertexIsNotEndpointOfParent flag
- if (secondaries) {
- size_t nseco = secondaries->size();
- for (size_t i = 0; i < nseco; i++) {
-
- // set current seco track ptr
- G4Track* seco = (*secondaries)[i];
-
- // setup new trk info for seco
- TrackInformation* secoTrkInfo = new TrackInformation(parTrkInfo);
- seco->SetUserInformation(secoTrkInfo);
-
- // set vtxnotendp flag in trk info
- if (parTrj) {
-
- // get parent's endpoint
- G4ThreeVector parEndpoint = getTrajectoryEndpoint(parTrj);
-
- // check if direct parent (flag always on if not) OR endp not near
- if (secoTrkInfo->getOriginalTrackID() != aTrack->GetTrackID()
- || vertexIsNotEndpointOfParent(parEndpoint, seco)) {
-
- // set vertexIsNotEndpointOfParent flag in trk info
- secoTrkInfo->setVertexIsNotEndpointOfParent(true);
- }
- }
- // vtxnotendp always set, because parent is not stored
- else {
- secoTrkInfo->setVertexIsNotEndpointOfParent(true);
- }
- }
-
-#ifdef SLIC_LOG
- log() << LOG::verbose << "track info passed to <" << nseco << "> secondaries" << LOG::endl;
-#endif
- }
-}
-
-bool TrajectoryManager::needsManualTrajectorySetup(const G4Track* aTrack) {
- // store != ON and has a corresponding tracker hit
- return (!m_trackingManager->GetStoreTrajectory())
- && (TrackInformation::getTrackInformation(aTrack)->hasTrackerHit());
-}
-
-void TrajectoryManager::stepping(const G4Step* aStep) {
- const G4Track* theTrack = aStep->GetTrack();
-
- // step has a live track
- if (isAlive(aStep)) {
-
- TrackInformation* trkInfo = TrackInformation::getTrackInformation(theTrack);
-
- UserRegionInformation* preReg = UserRegionInformation::getRegionInformation(
- aStep->GetPreStepPoint());
- UserRegionInformation* postReg = UserRegionInformation::getRegionInformation(
- aStep->GetPostStepPoint());
-
- // set tracking status based on region info
- if (preReg && postReg) {
-
- // tracking to non-tracking
- if (preReg->getStoreSecondaries() && (!postReg->getStoreSecondaries())) {
- trkInfo->setTrackingStatus(TrackInformation::eInNontrackingRegion);
- }
- // non-tracking to tracking
- else if (!preReg->getStoreSecondaries() && postReg->getStoreSecondaries()) {
- TrackInformation::ETrackingStatus origStatus = trkInfo->getOriginalTrackingStatus();
-
- // if looped like tracking -> nontracking -> tracking without dying
- if (!(origStatus == TrackInformation::eInTrackingRegion)) {
- trkInfo->setBackscatter(true);
- }
-
- // set tracking status to tracking
- trkInfo->setTrackingStatus(TrackInformation::eInTrackingRegion);
- }
- }
- }
-}
-
-void TrajectoryManager::setTrajectoryFinalStatus(const G4Track* aTrack) {
- // get current traj
- Trajectory* trj = getCurrentTrajectory();
-
- // must have trajectory to set status
- if (trj) {
-
- G4StepPoint* postStep = aTrack->GetStep()->GetPostStepPoint();
-
- // get track info
- TrackInformation* trkInfo = TrackInformation::getTrackInformation(aTrack);
-
- // has endpoint daughters
- G4int nSecoAtEnd = getNumberOfSecondariesAtEnd(aTrack);
- if (nSecoAtEnd > 0) {
- trj->setHasEndpointDaughters(true);
- }
-
- // default final status to unset
- Trajectory::EFinalStatus finalStatus = Trajectory::eUnset;
-
- // backscatter from track info
- trj->setBackscatter(trkInfo->getBackscatter());
-
- // vertexIsNotEndpointOfParent
- if (trkInfo->getVertexIsNotEndpointOfParent()) {
- trj->setVertexIsNotEndpointOfParent(true);
- }
-
- // no decay -> no secondaries at end
- if (nSecoAtEnd <= 0) {
- // left detector
- if (!postStep->GetMaterial()) {
- finalStatus = Trajectory::eLeftDetector;
- }
-
- // stopped
- else if (aTrack->GetKineticEnergy() <= 0.) {
- finalStatus = Trajectory::eStopped;
- }
- // unflagged loopers killed by G4
-#ifdef SLIC_DEBUG
- else {
- log() << LOG::debug << "Unknown decay case (probably a looper)." << LOG::done;
- }
-#endif
- }
- // decayed
- else {
-
- // decayed in tracker
- if (trkInfo->getTrackingStatus() == TrackInformation::eInTrackingRegion) {
- finalStatus = Trajectory::eDecayedInTracker;
- }
- // process of elimination: decayed in calorimeter
- else {
- finalStatus = Trajectory::eDecayedInCalorimeter;
- }
- }
-
- // set final status in trajectory
- trj->setFinalStatus(finalStatus);
- }
-}
-
-void TrajectoryManager::printTrack(const G4Track* trk, bool isSecondary) {
- log() << LOG::endl;
- log() << "track ptr <" << trk << ">" << LOG::endl;
- log() << "trackID <" << trk->GetTrackID() << ">" << LOG::endl;
- log() << "parent trackID <" << trk->GetParentID() << ">" << LOG::endl;
- log() << "PDG <" << trk->GetDefinition()->GetPDGEncoding() << ">" << LOG::endl;
- log() << "position <" << trk->GetPosition() << ">" << LOG::endl;
-
- if (!isSecondary) {
- if (trk->GetMaterial()) {
- log() << "material <" << trk->GetMaterial()->GetName() << ">" << LOG::endl;
- } else {
- log() << "NO material" << LOG::endl;
- }
-
- if (trk->GetVolume()) {
- log() << "volume <" << trk->GetVolume()->GetName() << ">" << LOG::endl;
- } else {
- log() << "NO volume" << LOG::endl;
- }
- }
-
- if (!isSecondary) {
- log() << "numSecoAtEnd <" << getNumberOfSecondariesAtEnd(trk) << ">" << LOG::endl;
- }
-
- log() << "kE <" << trk->GetKineticEnergy() << ">" << LOG::endl;
- log() << "mom <" << trk->GetMomentum() << ">" << LOG::endl;
-
- printStatusFlags(trk);
-
- log() << LOG::endl;
-}
-
-void TrajectoryManager::printPhysVolDbg(G4StepPoint* stepPnt) {
- if (stepPnt) {
- G4VPhysicalVolume* pv = stepPnt->GetPhysicalVolume();
-
- if (pv) {
- log() << LOG::debug << "PV name <" << pv->GetName() << ">" << LOG::endl;
- } else {
- log() << LOG::debug << "PV is NULL!" << LOG::endl;
- }
- }
-}
-
-G4int TrajectoryManager::getNumberOfSecondariesAtEnd(const G4Track* aTrack) {
- G4int numSecoAtEnd = 0;
- G4TrackVector* seco = m_trackingManager->GimmeSecondaries();
- G4ThreeVector endPnt = aTrack->GetPosition();
-
- if (seco) {
- size_t numSeco = seco->size();
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "numSeco <" << numSeco << ">" << LOG::endl;
-#endif
-
- for (size_t i = 0; i < numSeco; i++) {
- // if seco endpoint near track's endpoint
- if ((*seco)[i]->GetPosition().isNear(endPnt)) {
- ++numSecoAtEnd;
- }
- }
- }
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "found numSecoAtEnd <" << numSecoAtEnd << ">" << LOG::endl;
-#endif
-
- return numSecoAtEnd;
-}
-
-void TrajectoryManager::fillTrajectoryMap(const G4Event* anEvent, bool clearIt) {
-#ifdef SLIC_LOG
- log() << LOG::verbose << "fillTrackIDToTrajectoryMap()" << LOG::endl;
-#endif
- if (clearIt) {
- clearTrackIDToTrajectoryMap();
- }
-
- G4TrajectoryContainer* trajCont = anEvent->GetTrajectoryContainer();
-
- if (trajCont) {
- G4int n_traj = trajCont->entries();
-
- Trajectory* trj = 0;
- for (int i = 0; i < n_traj; i++) {
- trj = static_cast<Trajectory*>((*trajCont)[i]);
- addTrackIDToTrajectoryLink(trj);
- }
- } else {
- log() << LOG::error << "Cannot fill map without a trajectory container!" << LOG::endl;
- }
-}
-
-void TrajectoryManager::addTrackIDToTrajectoryLink(Trajectory* trj) {
- if (trj) {
- m_trackIDToTrajectory[trj->GetTrackID()] = trj;
- } else {
- log() << LOG::error << "Trajectory is null!" << LOG::endl;
- }
-}
-
-void TrajectoryManager::printSecondaries() {
- log() << LOG::debug << "secondaries for trackID <" << m_trackingManager->GetTrack()->GetTrackID() << ">"
- << LOG::done;
-
- G4TrackVector* secondaries = m_trackingManager->GimmeSecondaries();
- size_t nseco = secondaries->size();
- for (size_t i = 0; i < nseco; i++) {
-
- // set current seco track ptr
- G4Track* seco = (*secondaries)[i];
- log() << LOG::debug << "secondary track <" << i << ">" << LOG::endl;
- printTrack(seco, true);
- }
-}
-
-// find trajectory by track ID
-Trajectory* TrajectoryManager::findTrajectory(G4int trkID) {
-#ifdef SLIC_LOG
- log() << LOG::debug << "findTrajectory() - trkID <" << trkID << ">" << LOG::endl;
-#endif
-
- Trajectory* trj = 0;
- TrackIDToTrajectoryMap::iterator iter;
-
- iter = m_trackIDToTrajectory.find(trkID);
-
- if (iter != m_trackIDToTrajectory.end()) {
- trj = iter->second;
- }
-
-#ifdef SLIC_LOG
- if ( trj ) {
- log() << LOG::debug << "found trajectory" << LOG::endl;
- }
- else {
- log() << LOG::debug << "no trajectory found" << LOG::endl;
- }
-#endif
-
- return trj;
-}
-}
SVNspam 0.1