Commit in projects/slic/trunk on MAIN | |||
CMakeLists.txt | +22 | -23 | 3186 -> 3187 |
cmake/FindHEPPDT.cmake | +2 | -2 | 3186 -> 3187 |
include/CommandLineOption.hh | +2 | -2 | 3186 -> 3187 |
/CommandLineProcessor.hh | +4 | -4 | 3186 -> 3187 |
/CommandQueue.hh | +2 | -2 | 3186 -> 3187 |
/EventDebugger.hh | -129 | 3186 removed | |
/EventDebuggerMessenger.hh | -60 | 3186 removed | |
/EventSource.hh | +1 | -1 | 3186 -> 3187 |
/EventSourceManager.hh | +56 | -41 | 3186 -> 3187 |
/EventSourceWithInputFile.hh | +1 | -1 | 3186 -> 3187 |
/GPSEventSource.hh | +1 | -1 | 3186 -> 3187 |
/Geant4VersionInfo.hh.in | +30 | -14 | 3186 -> 3187 |
/GeneratorMessenger.hh | +1 | -1 | 3186 -> 3187 |
/LcioEventSource.hh | -92 | 3186 removed | |
/LcioHitsCollectionBuilder.hh | -3 | 3186 -> 3187 | |
/LcioManager.hh | +15 | -56 | 3186 -> 3187 |
/LcioMcpFactory.hh | -175 | 3186 removed | |
/LcioMcpFilter.hh | -69 | 3186 removed | |
/LcioMcpManager.hh | -230 | 3186 removed | |
/LcioMcpMaps.hh | -170 | 3186 removed | |
/LcioMcpMessenger.hh | -59 | 3186 removed | |
/LcioMcpPrinter.hh | -165 | 3186 removed | |
/LcioMcpStatusSetter.hh | -78 | 3186 removed | |
/LcioMcpUtil.hh | -54 | 3186 removed | |
/LcioMessenger.hh | +2 | -1 | 3186 -> 3187 |
/LcioPrimaryGenerator.hh | -76 | 3186 removed | |
/Module.hh | -14 | 3186 -> 3187 | |
/PackageInfo.hh.in | +19 | -6 | 3186 -> 3187 |
/ParticleGunEventSource.hh | +1 | -1 | 3186 -> 3187 |
/PrimaryGeneratorAction.hh | +7 | -17 | 3186 -> 3187 |
/RunManager.hh | -2 | 3186 -> 3187 | |
/SlicApplication.hh | +13 | -13 | 3186 -> 3187 |
/StdHepEventSource.hh | +4 | -12 | 3186 -> 3187 |
/StdHepLoader.hh | -98 | 3186 removed | |
/StdHepToLcioConvertor.hh | -195 | 3186 removed | |
/SteppingAction.hh | +7 | 3186 -> 3187 | |
/TrackingAction.hh | +7 | -9 | 3186 -> 3187 |
/Trajectory.hh | -337 | 3186 removed | |
/TrajectoryManager.hh | -405 | 3186 removed | |
macros/gps.mac | +1 | 3186 -> 3187 | |
src/EventAction.cc | +21 | -33 | 3186 -> 3187 |
/EventDebugger.cc | -113 | 3186 removed | |
/EventDebuggerMessenger.cc | -68 | 3186 removed | |
/EventSource.cc | +1 | -1 | 3186 -> 3187 |
/EventSourceManager.cc | +34 | -34 | 3186 -> 3187 |
/G4MagIntegratorStepperFactory.cc | +1 | 3186 -> 3187 | |
/GPSEventSource.cc | +1 | -1 | 3186 -> 3187 |
/GeneratorMessenger.cc | +5 | -4 | 3186 -> 3187 |
/LCExtendedParticles.cc | +3 | -7 | 3186 -> 3187 |
/LCOpticalPhysics.cc | +3 | -4 | 3186 -> 3187 |
/LCSUSYPhysics.cc | +75 | added 3187 | |
/LcioEventSource.cc | -93 | 3186 removed | |
/LcioFileNamer.cc | +5 | -2 | 3186 -> 3187 |
/LcioHitsCollectionBuilder.cc | +8 | -24 | 3186 -> 3187 |
/LcioManager.cc | +53 | -97 | 3186 -> 3187 |
/LcioMcpFactory.cc | -659 | 3186 removed | |
/LcioMcpFilter.cc | -173 | 3186 removed | |
/LcioMcpManager.cc | -264 | 3186 removed | |
/LcioMcpMaps.cc | -185 | 3186 removed | |
/LcioMcpMessenger.cc | -83 | 3186 removed | |
/LcioMcpPrinter.cc | -338 | 3186 removed | |
/LcioMcpStatusSetter.cc | -99 | 3186 removed | |
/LcioMcpUtil.cc | -31 | 3186 removed | |
/LcioMessenger.cc | +16 | -9 | 3186 -> 3187 |
/LcioPrimaryGenerator.cc | -261 | 3186 removed | |
/ParticleGunEventSource.cc | +1 | -1 | 3186 -> 3187 |
/PrimaryGeneratorAction.cc | +25 | -155 | 3186 -> 3187 |
/RunAction.cc | +6 | -2 | 3186 -> 3187 |
/RunManager.cc | +7 | -5 | 3186 -> 3187 |
/SlicApplication.cc | +8 | -63 | 3186 -> 3187 |
/StdHepEventSource.cc | +25 | -44 | 3186 -> 3187 |
/StdHepLoader.cc | -87 | 3186 removed | |
/StdHepToLcioConvertor.cc | -541 | 3186 removed | |
/SteppingAction.cc | +51 | -7 | 3186 -> 3187 |
/TrackingAction.cc | +80 | -17 | 3186 -> 3187 |
/Trajectory.cc | -127 | 3186 removed | |
/TrajectoryManager.cc | -584 | 3186 removed | |
+627 | -6834 |
Merge in v04-00-00-pre dev branch to 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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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/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; -} -}
Use REPLY-ALL to reply to list
To unsubscribe from the LCDET-SVN list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCDET-SVN&A=1