10 added + 22 removed + 20 modified, total 52 files
projects/slic/branches/v00-00-01_SLIC-218-dev
--- projects/slic/branches/v00-00-01_SLIC-218-dev/CMakeLists.txt 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/CMakeLists.txt 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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 )
@@ -76,17 +77,12 @@
# 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}" )
ENDIF()
-
-#EXECUTE_PROCESS( COMMAND grep "define G4VERSION_NUMBER" ${Geant4_DIR}/source/global/management/include/G4Version.hh | sed 's/#define G4VERSION_NUMBER\ \ //g' OUTPUT_VARIABLE GEANT4_VERSION )
-#EXECUTE_PROCESS( COMMAND grep "define G4VERSION_NUMBER" ${Geant4_DIR}/source/global/management/include/G4Version.hh OUTPUT_VARIABLE GEANT4_VERSION )
-#MESSAGE( STATUS "GEANT4_VERSION=${GEANT4_VERSION}" )
-
# 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 )
@@ -100,6 +96,9 @@
ADD_DEFINITIONS( -DSLIC_LOG )
ENDIF()
+# debug flags
+SET ( CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS} -DDEBUG=1 -pg -ggdb" )
+
# local includes
SET( SLIC_INCLUDE_DIR ${PROJECT_SOURCE_DIR}/include )
INCLUDE_DIRECTORIES( ${SLIC_INCLUDE_DIR} )
@@ -155,8 +154,7 @@
ENDIF()
-
-#get_cmake_property(_variableNames VARIABLES)
-#foreach (_variableName ${_variableNames})
-# message(STATUS "${_variableName}=${${_variableName}}")
-#endforeach()
+get_cmake_property(_variableNames VARIABLES)
+foreach (_variableName ${_variableNames})
+ message(STATUS "${_variableName}=${${_variableName}}")
+endforeach()
projects/slic/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/EventSourceManager.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/EventSourceManager.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -291,44 +291,44 @@
private:
// generator messenger
- GeneratorMessenger* m_messenger;
+ GeneratorMessenger* m_messenger;
- // current event source
- EventSource* m_currentEventSource;
+ // current event source
+ EventSource* m_currentEventSource;
- // ParticleGun source. Always enabled.
- EventSource* m_particleGunSource;
+ // ParticleGun source. Always enabled.
+ EventSource* m_particleGunSource;
- // filename
- std::string m_filename;
- bool m_fileIsSet;
- bool m_newFilename;
+ // filename
+ std::string m_filename;
+ bool m_fileIsSet;
+ bool m_newFilename;
- // num events generated on current generator
- int m_ngen;
+ // num events generated on current generator
+ int m_ngen;
- // num events to skip
- unsigned int m_nskip;
+ // num events to skip
+ unsigned int m_nskip;
- // need to setup a new source new run?
- bool m_newSource;
+ // need to setup a new source new run?
+ bool m_newSource;
- // gen strings
- static std::string m_stdhepStr;
- static std::string m_lcioStr;
- static std::string m_gpsStr;
- static std::string m_gunStr;
- static std::string m_unknownStr;
+ // 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;
- // source type
- ESourceType m_sourceType;
+ // source type
+ ESourceType m_sourceType;
- // Lorentz transformation angle which will be applied to event.
- G4double m_lorentzTransformationAngle;
+ // 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;
- };
- }
+ // Parameter which will be used to smear vertex Z position.
+ G4double m_zSmearingParam;
+};
+}
#endif
projects/slic/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioEventSource.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioEventSource.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioHitsCollectionBuilder.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioHitsCollectionBuilder.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -193,9 +193,6 @@
// store momentum
bool m_storeMomentum;
-
- // ptr to Lcio Mcp manager for lkp of Mcp data
- LcioMcpManager* m_mcpManager;
};
}
projects/slic/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioManager.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioManager.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -9,8 +9,6 @@
#include "lcdd/hits/TrackerHit.hh"
// slic
-#include "LcioPrimaryGenerator.hh"
-#include "LcioMcpManager.hh"
#include "LcioHitsCollectionBuilder.hh"
#include "LcioFileNamer.hh"
#include "Trajectory.hh"
@@ -233,22 +231,6 @@
}
/**
- * 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.
*/
@@ -327,9 +309,6 @@
private:
- // final Mcp collection with all info from McpInitial, Trajectories, Primaries
- IMPL::LCCollectionVec* m_McpFinalColl;
-
// writer
IO::LCWriter* m_writer;
@@ -340,11 +319,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;
projects/slic/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpFactory.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpFactory.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpFilter.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpFilter.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpManager.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpManager.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpMaps.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpMaps.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpMessenger.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpMessenger.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpPrinter.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpPrinter.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpStatusSetter.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpStatusSetter.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpUtil.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioMcpUtil.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioPrimaryGenerator.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/LcioPrimaryGenerator.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/MCParticleManager.hh (rev 0)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/MCParticleManager.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -0,0 +1,65 @@
+#ifndef slic_MCParticleManager_hh_
+#define slic_MCParticleManager_hh_ 1
+
+// SLIC
+#include "Singleton.hh"
+
+// LCIO
+#include "IMPL/LCCollectionVec.h"
+#include "IMPL/MCParticleImpl.h"
+#include "EVENT/MCParticle.h"
+
+// STL
+#include <map>
+
+using std::map;
+using EVENT::MCParticle;
+using IMPL::MCParticleImpl;
+using IMPL::LCCollectionVec;
+
+namespace slic {
+
+typedef map<G4int, MCParticle*> MCParticle2TrackIDMap_type;
+
+class MCParticleManager : public Singleton<MCParticleManager> {
+
+public:
+
+ MCParticleManager()
+ : _mcpVec(0) {
+ }
+
+ virtual ~MCParticleManager() {}
+
+ LCCollectionVec* getMCParticleCollection() {
+ return _mcpVec;
+ }
+
+ void setMCParticleCollection(LCCollectionVec* mcpVec) {
+ _mcpVec = mcpVec;
+ }
+
+ void addMCParticleTrackID(G4int trackId, MCParticle* particle) {
+ if (trackId <= 0) {
+ return;
+ }
+ _mcpMap[trackId] = particle;
+ }
+
+ MCParticle* findMCParticle(G4int trackId) {
+ return _mcpMap[trackId];
+ }
+
+ void clearMap() {
+ _mcpMap.clear();
+ }
+
+private:
+
+ MCParticle2TrackIDMap_type _mcpMap;
+ LCCollectionVec* _mcpVec;
+};
+
+}
+
+#endif
projects/slic/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/PrimaryGeneratorAction.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/PrimaryGeneratorAction.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -46,14 +46,14 @@
* is retrieved from the singleton EventSourceManager.
* @param[in] anEvent The current Geant4 event.
*/
- void applyLorentzTransformation(G4Event *anEvent);
+ //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);
+ //void zsmear(G4Event* anEvent);
private:
projects/slic/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/StdHepEventSource.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/StdHepEventSource.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -5,8 +5,8 @@
// slic
#include "StdHepLoader.hh"
-#include "StdHepToLcioConvertor.hh"
#include "EventSourceWithInputFile.hh"
+#include "StdHepGenerator.hh"
namespace slic {
/**
@@ -72,7 +72,7 @@
private:
StdHepLoader* m_loader;
- StdHepToLcioConvertor* m_convertor;
+ StdHepGenerator* m_eventGenerator;
};
}
projects/slic/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/StdHepGenerator.hh (rev 0)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/StdHepGenerator.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -0,0 +1,57 @@
+#ifndef slic_StdhepGenerator_hh_
+#define slic_StdhepGenerator_hh_ 1
+
+// LCIO
+#include "UTIL/LCStdHepRdr.h"
+#include "EVENT/MCParticle.h"
+
+// Geant4
+#include "G4Event.hh"
+#include "G4PrimaryParticle.hh"
+
+// STL
+#include <set>
+
+using UTIL::LCStdHepRdr;
+using EVENT::MCParticle;
+using EVENT::LCCollection;
+
+namespace slic {
+
+/**
+ * @class StdHepGenerator
+ * @brief Converts from StdHep input to a G4Event
+ * @note Borrows heavily from Mokka's HepLCIOInterfaceNew class.
+ */
+class StdHepGenerator {
+
+public:
+
+ typedef std::map<MCParticle*, G4PrimaryParticle*> LCIO2Geant4Map;
+ typedef LCIO2Geant4Map::iterator LCIO2Geant4MapIterator;
+
+public:
+
+ StdHepGenerator(G4String eventFile);
+
+ virtual ~StdHepGenerator();
+
+ void generateEvent(G4Event*);
+
+ static LCIO2Geant4Map* getParticleMap() {
+ return &_particleMap;
+ }
+
+private:
+
+ std::set<G4PrimaryParticle*> getRelevantParticles(MCParticle* p);
+
+private:
+ static LCIO2Geant4Map _particleMap;
+ static std::set<MCParticle*> _visitedParticles;
+ LCStdHepRdr* _reader;
+};
+
+}
+
+#endif
projects/slic/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/StdHepToLcioConvertor.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/StdHepToLcioConvertor.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/TrackManager.hh (rev 0)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/TrackManager.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -0,0 +1,68 @@
+#ifndef slic_TrackManager_hh_
+#define slic_TrackManager_hh_ 1
+
+// SLIC
+#include "Module.hh"
+#include "TrackSummary.hh"
+#include "Singleton.hh"
+
+// LCIO
+#include "EVENT/LCEvent.h"
+
+// Geant4
+#include "G4Event.hh"
+
+using EVENT::LCEvent;
+
+namespace slic {
+
+/**
+ * @class TrackManager
+ * @brief Manages access to stored track information used for MCParticle output
+ */
+class TrackManager : public Singleton<TrackManager>, public Module {
+
+public:
+
+ TrackManager() : Module("TrackManager") {
+ _trackSummaries = new TrackSummaryVector;
+ }
+
+ virtual ~TrackManager() {
+ delete _trackSummaries;
+ }
+
+public:
+
+ TrackSummaryVector* getTrackSummaries() {
+ return _trackSummaries;
+ }
+
+ void resetTrackSummaries() {
+ _trackSummaries->clear();
+ }
+
+ void saveTracks(const G4Event*, LCEvent*);
+
+ void clearTrackSummaries() {
+ for (TrackSummaryVector::iterator it = _trackSummaries->begin();
+ it != _trackSummaries->end();
+ it++) {
+ delete *it;
+ }
+ _trackSummaries->clear();
+ }
+
+private:
+
+ TrackSummary* findFirstSavedAncestor(TrackSummary* trackSummary);
+
+private:
+
+ TrackSummaryVector* _trackSummaries;
+
+};
+
+}
+
+#endif
projects/slic/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/TrackSummary.hh (rev 0)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/TrackSummary.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -0,0 +1,169 @@
+#ifndef slic_TrackSummary_hh
+#define slic_TrackSummary_hh 1
+
+// LCIO
+#include "lcio.h"
+#include "IMPL/MCParticleImpl.h"
+
+// Geant4
+#include "globals.hh"
+#include "G4ThreeVector.hh"
+#include "G4Allocator.hh"
+
+// STL
+#include <vector>
+#include <bitset>
+
+using IMPL::MCParticleImpl;
+
+class G4Track;
+
+namespace slic {
+
+class TrackSummary;
+typedef std::vector<TrackSummary*> TrackSummaryVector;
+
+class TrackSummary {
+
+public:
+
+ // Constructor/destructor
+ TrackSummary(const G4Track*, G4bool, G4bool backScattering = false);
+
+ virtual ~TrackSummary() {
+ }
+
+ // new and delete for allocator usage
+ inline void *operator new(size_t);
+ inline void operator delete(void *);
+
+ void update(const G4Track*);
+
+ //
+ // Get/Set functions
+ //
+
+ G4double getCharge() const {
+ return _charge;
+ }
+
+ G4ThreeVector getEndPoint() const {
+ return _endPoint;
+ }
+
+ G4double getEnergy() const {
+ return sqrt(_momentum[0] * _momentum[0] + _momentum[1] * _momentum[1] + _momentum[2] * _momentum[2] + _mass * _mass);
+ }
+
+ G4double getMass() const {
+ return _mass;
+ }
+
+ G4int getHepEvtStatus() const {
+ return _hepEvtStatus;
+ }
+
+ G4int getSimulatorStatus() const {
+ return _simstatus.to_ulong();
+ }
+
+ G4ThreeVector getMomentum() const {
+ return _momentum;
+ }
+
+ TrackSummary* getParent() const {
+ return _parent;
+ }
+
+ G4int getPDG() const {
+ return _PDG;
+ }
+
+ const G4ThreeVector& getVertex() const {
+ return _vertex;
+ }
+
+ G4int getTrackID() const {
+ return _trackID;
+ }
+
+ G4int getParentID() const {
+ return _parentTrackID;
+ }
+
+ G4double getTrackLength() const {
+ return _trackLength;
+ }
+
+ G4bool getToBeSaved() const {
+ return _save;
+ }
+
+ void setToBeSaved() {
+ _save = true;
+ }
+
+ G4bool getBackScattering() const {
+ return _backScattering;
+ }
+
+ void setBackScattering() {
+ _backScattering = true;
+ }
+
+ void setParentToBeSaved();
+
+ MCParticleImpl* getMCParticle() {
+ return _theMCParticle;
+ }
+
+ void setMCParticle(MCParticleImpl* mcp) {
+ _theMCParticle = mcp;
+ }
+
+ void buildMCParticle();
+
+ TrackSummary* findParent() const;
+
+private:
+
+// Member data
+
+ G4double _charge;
+ G4ThreeVector _endPoint;
+ G4double _energy;
+ G4double _mass;
+ G4int _hepEvtStatus;
+ std::bitset<32> _simstatus;
+ G4ThreeVector _momentum;
+ TrackSummary* _parent;
+ G4int _PDG;
+ G4ThreeVector _vertex;
+
+ G4bool _save;
+
+ G4bool _backScattering;
+ G4int _trackID;
+ G4int _parentTrackID;
+ G4double _trackLength;
+
+ MCParticleImpl* _theMCParticle;
+ G4bool _mcParticleIsUpToDate;
+
+};
+
+extern G4Allocator<TrackSummary> TrackSummaryAllocator;
+
+inline void* TrackSummary::operator new(size_t) {
+ void *aTrackSummary;
+ aTrackSummary = (void*) TrackSummaryAllocator.MallocSingle();
+ return aTrackSummary;
+}
+
+inline void TrackSummary::operator delete(void* aTrackSummary) {
+ TrackSummaryAllocator.FreeSingle((TrackSummary*) aTrackSummary);
+}
+
+}
+
+#endif
projects/slic/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/TrackUtil.hh (rev 0)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/TrackUtil.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -0,0 +1,61 @@
+#ifndef slic_TrackUtil_hh_
+#define slic_TrackUtil_hh_ 1
+
+// LCDD
+#include "lcdd/geant4/UserRegionInformation.hh"
+
+// SLIC
+#include "TrackManager.hh"
+#include "UserTrackInformation.hh"
+
+// Geant4
+#include "G4TransportationManager.hh"
+
+namespace slic {
+
+/**
+ * @class TrackUtil
+ * @brief Miscellaneous helper methods for G4Track and related classes
+ */
+class TrackUtil {
+
+public:
+
+ /**
+ * Get the UserRegionInformation of this G4Track based on its start position.
+ * @param track The G4Track.
+ * @return The UserRegionInformation for this G4Track.
+ */
+ static UserRegionInformation* getRegionInformation(const G4Track* track) {
+ G4VPhysicalVolume* volume = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->LocateGlobalPointAndSetup(
+ track->GetPosition());
+ G4Region* region = volume->GetLogicalVolume()->GetRegion();
+ UserRegionInformation* regionInfo = dynamic_cast<UserRegionInformation*>(region->GetUserInformation());
+ return regionInfo;
+ }
+
+ /**
+ * Setup the UserTrackInformation for the G4Track, which also requires first
+ * creating a TrackSummary containing all the track information.
+ * @return The UserTrackInformation that was created for the track.
+ */
+ static UserTrackInformation* setupUserTrackInformation(const G4Track* track, G4bool save) {
+ TrackSummary* trackSummary = new TrackSummary(track, save);
+ TrackManager::instance()->getTrackSummaries()->push_back(trackSummary);
+ UserTrackInformation* trackInfo = new UserTrackInformation(trackSummary);
+ ((G4Track*) track)->SetUserInformation(trackInfo);
+ return trackInfo;
+ }
+
+ /**
+ * Get the UserTrackInformation for the track.
+ * @return The UserTrackInformation for the track, which may be NULL.
+ */
+ static UserTrackInformation* getUserTrackInformation(const G4Track* track) {
+ return dynamic_cast<UserTrackInformation*>(track->GetUserInformation());
+ }
+};
+
+};
+
+#endif
projects/slic/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/TrackingAction.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/TrackingAction.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/TrajectoryManager.hh 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/TrajectoryManager.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -356,6 +356,14 @@
}
/**
+ * Return true if the G4Track is currently stepping within the world volume.
+ * @return True if G4Track is within the world volume; false if in a sub-volume.
+ */
+ static bool inWorldVolume(const G4Track* track) {
+ return track->GetStep()->GetPreStepPoint()->GetPhysicalVolume()->GetMotherLogical() == 0;
+ }
+
+ /**
* Print the status flags of a G4Track.
* @param[in] track The G4Track.
*/
projects/slic/branches/v00-00-01_SLIC-218-dev/include
--- projects/slic/branches/v00-00-01_SLIC-218-dev/include/UserTrackInformation.hh (rev 0)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/include/UserTrackInformation.hh 2014-01-09 20:37:54 UTC (rev 2864)
@@ -0,0 +1,59 @@
+#ifndef slic_UserTrackInformation_hh_
+#define slic_UserTrackInformation_hh_ 1
+
+// SLIC
+#include "TrackSummary.hh"
+
+// Geant4
+#include "globals.hh"
+#include "G4Allocator.hh"
+#include "G4VUserTrackInformation.hh"
+
+namespace slic {
+
+/**
+ * @class UserTrackInformation
+ * @brief The implementation of the G4VUserTrackInformation class for storing track information within the event
+ */
+class UserTrackInformation: public G4VUserTrackInformation {
+
+public:
+
+ UserTrackInformation(TrackSummary * aTrackSummary) {
+ _trackSummary = aTrackSummary;
+ }
+
+ virtual ~UserTrackInformation() {
+ }
+
+ inline void *operator new(size_t);
+
+ inline void operator delete(void *trackInfo);
+
+ TrackSummary* getTrackSummary() const {
+ return _trackSummary;
+ }
+
+ void Print() const {
+ }
+
+private:
+ TrackSummary* _trackSummary;
+};
+
+extern G4Allocator<UserTrackInformation> TrackInformationAllocator;
+
+inline void* UserTrackInformation::operator new(size_t) {
+ void* trackInfo;
+ trackInfo = (void*) TrackInformationAllocator.MallocSingle();
+ return trackInfo;
+}
+
+inline void UserTrackInformation::operator delete(void *trackInfo) {
+ TrackInformationAllocator.FreeSingle((UserTrackInformation*) trackInfo);
+}
+
+}
+
+#endif
+
projects/slic/branches/v00-00-01_SLIC-218-dev/macros
--- projects/slic/branches/v00-00-01_SLIC-218-dev/macros/gps.mac 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/macros/gps.mac 2014-01-09 20:37:54 UTC (rev 2864)
@@ -1,3 +1,4 @@
+/generator/select gps
/gps/verbose 2
/gps/particle mu-
/gps/energy 10 GeV
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/EventAction.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/EventAction.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -5,12 +5,14 @@
#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 "EventDebugger.hh"
+#include "LcioManager.hh"
+#include "MCParticleManager.hh"
+#include "SlicApplication.hh"
+#include "TrackManager.hh"
+#include "TrajectoryManager.hh"
// Geant4
#include "G4RunManager.hh"
@@ -28,36 +30,32 @@
}
void EventAction::BeginOfEventAction(const G4Event *anEvent) {
- // check aborted
- if (!anEvent->IsAborted()) {
- // debugger
- //EventDebugger::instance()->beginEvent( anEvent );
+ /* Clear the MCParticle map from last event. */
+ MCParticleManager::instance()->clearMap();
- // trj mgr
- TrajectoryManager::instance()->beginEvent(anEvent);
-
- /* 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) {
+
+ /* Save track information to MCParticle collection in LCIO output event. */
+ TrackManager::instance()->saveTracks(anEvent, LcioManager::instance()->getCurrentLCEvent());
+
// check aborted
if (!SlicApplication::instance()->isAborting()) {
- // LcioManager's action
+ // Execute LcioManager's end of event action.
LcioManager::instance()->endEvent(anEvent);
}
- // event source (generator) action
+ /* Begin of event processing for the current event generator. */
EventSourceManager::instance()->endEvent(anEvent);
- // event timer
+ /* Clear the track information for this event. */
+ TrackManager::instance()->clearTrackSummaries();
+
+ /* Stop the event timer. */
stopEventTimer();
// debugger
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/EventSourceManager.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/EventSourceManager.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -1,15 +1,14 @@
// $Header: /nfs/slac/g/lcd/cvs/lcdroot/slic/src/EventSourceManager.cc,v 1.26 2013-06-26 01:57:02 jeremy Exp $
// SLIC
+#include "EventAction.hh"
#include "EventSourceManager.hh"
+#include "FileUtil.hh"
#include "GPSEventSource.hh"
#include "ParticleGunEventSource.hh"
+#include "RunManager.hh"
+#include "SlicApplication.hh"
#include "StdHepEventSource.hh"
-#include "LcioEventSource.hh"
-#include "FileUtil.hh"
-#include "EventAction.hh"
-#include "SlicApplication.hh"
-#include "RunManager.hh"
// LCDD
#include "lcdd/util/StringUtil.hh"
@@ -23,9 +22,9 @@
std::string EventSourceManager::m_unknownStr = 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.) {
// messenger with generator command macro bindings
m_messenger = new GeneratorMessenger();
@@ -175,9 +174,10 @@
EventSource* EventSourceManager::createEventSource(ESourceType st) {
EventSource* src = 0;
- if (st == eLCIO) {
+ /*if (st == eLCIO) {
src = new LcioEventSource(getFilename());
- } else if (st == eStdHep) {
+ } else*/
+ if (st == eStdHep) {
src = new StdHepEventSource(getFilename());
} else if (st == eGPS) {
src = new GPSEventSource();
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioEventSource.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioEventSource.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioHitsCollectionBuilder.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioHitsCollectionBuilder.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -2,8 +2,8 @@
// SLIC
#include "LcioHitsCollectionBuilder.hh"
-#include "LcioMcpManager.hh"
#include "HitsCollectionUtil.hh"
+#include "MCParticleManager.hh"
#include "SlicApplication.hh"
// LCDD
@@ -13,18 +13,18 @@
#include "G4SDManager.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>
@@ -35,8 +35,6 @@
LcioHitsCollectionBuilder::LcioHitsCollectionBuilder() :
Module("LcioHitsCollectionBuilder"), m_storeMomentum(0) {
- // set local LcioMcpManager ptr
- m_mcpManager = LcioMcpManager::instance();
// setup default coll flag for cal hits
setFlagDefaults();
@@ -287,20 +285,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);
}
@@ -331,37 +323,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/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioManager.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioManager.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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"
-// 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;
@@ -54,18 +53,11 @@
string LcioManager::m_defaultFileName = "outfile";
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();
@@ -376,20 +368,6 @@
void LcioManager::endEvent(const G4Event*) {
if (!RunManager::instance()->isRunAborted()) {
- // create LCEvent
- createLCEvent();
-
- // create Mcp coll in LcioMcpManager
- m_mcpManager->endEvent(G4EventManager::GetEventManager()->GetNonconstCurrentEvent());
-
- // create mcp collection from input event
- createFinalMcpCollection();
-
- // If selected, add the initial MCParticle collection to the event.
- if (LcioMcpManager::instance()->writeInitialMCParticleCollection()) {
- addInitialMCParticleCollection();
- }
-
// create HC in current LCEvent from current G4Event using builder
createHitsCollections();
@@ -425,16 +403,8 @@
// 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 +414,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/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpFactory.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpFactory.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -1,658 +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"
-
-// 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/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpFilter.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpFilter.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpManager.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpManager.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpMaps.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpMaps.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpMessenger.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpMessenger.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpPrinter.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpPrinter.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpStatusSetter.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpStatusSetter.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpUtil.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMcpUtil.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMessenger.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioMessenger.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -3,7 +3,6 @@
// SLIC
#include "LcioMessenger.hh"
#include "LcioManager.hh"
-#include "LcioMcpManager.hh"
#include "LcioFileNamer.hh"
#include "TrajectoryManager.hh"
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioPrimaryGenerator.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/LcioPrimaryGenerator.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -1,260 +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"
-
-// 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/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/PrimaryGeneratorAction.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/PrimaryGeneratorAction.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -5,8 +5,8 @@
#include "LcioManager.hh"
#include "RunManager.hh"
#include "SlicApplication.hh"
-#include "LcioMcpManager.hh"
-#include "LcioMcpFilter.hh"
+//#include "LcioMcpManager.hh"
+//#include "LcioMcpFilter.hh"
#include "EventSourceManager.hh"
// LCDD
@@ -30,6 +30,22 @@
void PrimaryGeneratorAction::GeneratePrimaries(G4Event *anEvent) {
+ /* Print begin of event message. */
+ printBeginEventMessage(anEvent);
+
+ /* Get the event source manager. */
+ EventSourceManager* mgr = EventSourceManager::instance();
+
+ /* Begin event hook of event source manager. */
+ mgr->beginEvent(anEvent);
+
+ /* Generate primary vertices using the current event generator. */
+ mgr->GeneratePrimaryVertex(anEvent);
+
+ // DEBUG: stop program here!!!
+ //assert(false);
+
+ /*
// Check if run needs to be aborted, in which case we immediately return without executing anything.
if (RunManager::instance()->isRunAborted()) {
G4RunManager::GetRunManager()->AbortRun();
@@ -37,22 +53,6 @@
return;
}
- EventSourceManager* mgr = EventSourceManager::instance();
-
- // 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);
-
- // 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();
-
- // Generator manager begin event hook.
- mgr->beginEvent(anEvent);
-
- // Generate vertices using the manager.
- mgr->GeneratePrimaryVertex(anEvent);
-
// Apply Lorentz transformation to event, used for non-zero beam crossing angle.
applyLorentzTransformation(anEvent);
@@ -63,6 +63,7 @@
if (mgr->isEOF()) {
SlicApplication::instance()->setAborting(true);
}
+ */
}
void PrimaryGeneratorAction::printBeginEventMessage(G4Event* anEvent) {
@@ -70,6 +71,9 @@
<< LOG::done;
}
+// TODO: The lorentz transform and z smearing functions below should be moved elsewhere!!!
+
+/*
void PrimaryGeneratorAction::applyLorentzTransformation(G4Event *evt) {
const G4double alpha = EventSourceManager::instance()->getLorentzTransformationAngle();
@@ -137,10 +141,12 @@
}
}
}
+*/
/**
* Apply a Z smear to the position of every G4PrimaryVertex in an event.
*/
+/*
void PrimaryGeneratorAction::zsmear(G4Event* anEvent) {
//G4cout << "zsmear" << G4endl;
@@ -187,4 +193,5 @@
}
}
}
+*/
}
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/RunAction.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/RunAction.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -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();
@@ -30,6 +33,9 @@
printBeginOfRunMessage( aRun );
#endif
+ // Reset the global TrackSummaryVector used for MCParticle output.
+ TrackManager::instance()->resetTrackSummaries();
+
// Execute LcioManager's beginRun action.
LcioManager::instance()->beginRun(aRun);
@@ -44,9 +50,14 @@
}
void RunAction::EndOfRunAction(const G4Run *run) {
+
// LcioManager's action
LcioManager::instance()->endRun(run);
+ std::cout << "TrackManager has "
+ << TrackManager::instance()->getTrackSummaries()->size()
+ << " TrackSummary objects at end of run." << std::endl;
+
// event source mgr
EventSourceManager::instance()->endRun(run);
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/StdHepEventSource.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/StdHepEventSource.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -1,23 +1,20 @@
// $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"
namespace slic {
StdHepEventSource::StdHepEventSource(const std::string& fname) :
- EventSourceWithInputFile("StdHepEventSource", fname), m_convertor(0) {
+ EventSourceWithInputFile("StdHepEventSource", fname),
+ m_eventGenerator(new StdHepGenerator(fname)) {
+
// create loader with internal reader
m_loader = new StdHepLoader();
-
- // create new convertor
- m_convertor = new StdHepToLcioConvertor();
}
void StdHepEventSource::generate(G4Event* anEvent) {
- LcioManager::instance()->getGenerator()->generatePrimaryVertexFromMcpCollection(
- LcioMcpManager::instance()->getInitialMcpCollection(), anEvent);
+ m_eventGenerator->generateEvent(anEvent);
}
// open the current file
@@ -60,23 +57,12 @@
EventSourceWithInputFile::beginRun(aRun);
// setup convertor
- assert( m_loader->getStdHepReader());
- m_convertor->setStdHepReader(m_loader->getStdHepReader());
+ assert(m_loader->getStdHepReader());
}
void StdHepEventSource::beginEvent(const G4Event* anEvent) {
// read an event
EventSourceWithInputFile::beginEvent(anEvent);
-
- // not EOF?
- if (!isEOF()) {
-
- // set convertor's coll
- m_convertor->setParticleCollection(LcioMcpManager::instance()->getInitialMcpCollection());
-
- // run convertor
- m_convertor->fillParticleCollection();
- }
}
lStdHep* StdHepEventSource::getStdHepReader() {
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/StdHepGenerator.cc (rev 0)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/StdHepGenerator.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -0,0 +1,206 @@
+#include "StdHepGenerator.hh"
+
+// SLIC
+#include "MCParticleManager.hh"
+
+// LCIO
+#include "EVENT/MCParticle.h"
+#include "IMPL/MCParticleImpl.h"
+
+using UTIL::LCStdHepRdr;
+using EVENT::MCParticle;
+using EVENT::LCCollection;
+using IMPL::MCParticleImpl;
+
+namespace slic {
+
+StdHepGenerator::LCIO2Geant4Map StdHepGenerator::_particleMap;
+std::set<MCParticle*> StdHepGenerator::_visitedParticles;
+
+StdHepGenerator::StdHepGenerator(G4String eventFile) {
+ _reader = new LCStdHepRdr(eventFile.data());
+}
+
+StdHepGenerator::~StdHepGenerator() {
+ delete _reader;
+}
+
+void StdHepGenerator::generateEvent(G4Event* event) {
+
+ /* Clear previous state before starting event generation. */
+ _particleMap.clear();
+ _visitedParticles.clear();
+
+ /* Get collection of MCParticles from the LCStdHepRdr. */
+ LCCollectionVec* mcpVec = _reader->readEvent();
+ std::vector<MCParticle*> collection(mcpVec->getNumberOfElements());
+ G4int n = collection.size();
+
+ /*
+ * Copy input MCParticles to a vector in order to work with them easily.
+ */
+ for(G4int i=0; i<n; i++) {
+ MCParticleImpl* mcp = dynamic_cast<MCParticleImpl*>(mcpVec->getElementAt(i));
+ collection[i] = mcp;
+ }
+
+ /*
+ * Process particles and add them to the primary vertex of the event.
+ */
+ G4ThreeVector particlePosition; // FIXME: add Z smear
+ G4double particleTime;
+ G4PrimaryVertex* vertex = new G4PrimaryVertex(particlePosition, particleTime);
+ for(size_t i=0; i < collection.size(); i++) {
+ MCParticle* mcp = collection[i];
+ if (mcp->getParents().size()==0) {
+ //std::cout << "construct tree of relevant G4PrimaryParticles for initial MCParticle " << i << " " << std::flush;
+ std::set<G4PrimaryParticle*> primaries = getRelevantParticles(mcp);
+ //std::cout << "... finished" << std::endl;
+ std::set<G4PrimaryParticle*>::iterator primaryIt;
+ for (primaryIt = primaries.begin(); primaryIt != primaries.end(); primaryIt++) {
+ vertex->SetPrimary(*primaryIt);
+ //std::cout << "G4PrimaryParticle ("<< (*primaryIt)->GetPDGcode() << ") added to G4PrimaryVertex." << std::endl;
+ }
+ }
+ }
+
+ /* Add the primary vertex to Geant4. */
+ event->AddPrimaryVertex(vertex);
+
+ /* Register the collection of MCParticles that was created. */
+ MCParticleManager::instance()->setMCParticleCollection(mcpVec);
+}
+
+std::set<G4PrimaryParticle*> StdHepGenerator::getRelevantParticles(MCParticle* p) {
+
+ //log each particle which has been called, to avoid double counting and increase efficiency
+ _visitedParticles.insert(p);
+ LCIO2Geant4Map::iterator mcpIT;
+ std::set<G4PrimaryParticle*> relevantParticlesSet; //holds all relevant decay particles of p
+
+ /*
+ * Case #1:
+ *
+ * If p is a stable particle: search for it in LCIO2Geant4Map
+ *
+ * If exists in map: get the G4PrimaryParticle version of p from LCIO2Geant4Map
+ * Else: create G4PrimaryParticle version write it to LCIO2Geant4Map
+ *
+ * Finally: insert the G4PrimaryParticle verion of p into the relevant particle list
+ * and return the list.
+ */
+ if (p->getGeneratorStatus() == 1) {
+ G4PrimaryParticle* g4p;
+ mcpIT = _particleMap.find(p);
+ if (mcpIT != _particleMap.end()) {
+ g4p = mcpIT->second;
+ } else {
+ G4int IDHEP = p->getPDG();
+ G4double PHEP1 = p->getMomentum()[0];
+ G4double PHEP2 = p->getMomentum()[1];
+ G4double PHEP3 = p->getMomentum()[2];
+ G4double PHEP5 = p->getMass();
+ // create G4PrimaryParticle object
+ g4p = new G4PrimaryParticle(IDHEP, PHEP1 * GeV, PHEP2 * GeV, PHEP3 * GeV);
+ g4p->SetMass(PHEP5 * GeV);
+ _particleMap[p] = g4p;
+ //std::cout << "*" << std::flush;
+ }
+ //std::cout << g4p->GetPDGcode() << std::flush;
+ relevantParticlesSet.insert(g4p);
+ return relevantParticlesSet;
+ }
+
+ /*
+ * CASE #2:
+ *
+ * If p is not stable: get first decay daughter and calculate the proper time of p.
+ *
+ * If proper time is not zero: particle is "relevant", since it carries vertex information.
+ *
+ * If p is already in LCIO2Geant4Map: take it.
+ * Else: create G4PrimaryParticle version of p and write it to LCIO2Geant4Map
+ *
+ * Collect all relevant particles of all daughters and setup "relevant mother-
+ * daughter-relations" between relevant decay particles and G4PrimaryParticle version of p
+ *
+ * Finally: insert only the G4PrimaryParticle version of p to the relevant particle list and return this "list".
+ * The added particle has now the correct pre-assigned decay products and (hopefully) the right lifetime.
+ */
+ else if (p->getDaughters().size() > 0) { //fg: require that there is at least one daughter - otherwise forget the particle
+ // calculate proper time
+ MCParticle* dp = p->getDaughters()[0];
+
+ double proper_time = fabs((dp->getTime() - p->getTime()) * p->getMass()) / p->getEnergy();
+ // fix by S.Morozov for real != 0
+ double aPrecision = dp->getTime() * p->getMass() / p->getEnergy();
+ double bPrecision = p->getTime() * p->getMass() / p->getEnergy();
+
+ double proper_time_Precision = pow(10, -DBL_DIG)*fmax(fabs(aPrecision),fabs(bPrecision));
+
+ bool isProperTimeZero = (proper_time <= proper_time_Precision);
+
+ // -- remove original --- if (proper_time != 0) {
+ if (isProperTimeZero == false) {
+
+ G4PrimaryParticle* g4p;
+ mcpIT = _particleMap.find(p);
+ if (mcpIT != _particleMap.end()) {
+ g4p = mcpIT->second;
+ } else {
+ G4int IDHEP = p->getPDG();
+ G4double PHEP1 = p->getMomentum()[0];
+ G4double PHEP2 = p->getMomentum()[1];
+ G4double PHEP3 = p->getMomentum()[2];
+ G4double PHEP5 = p->getMass();
+ // create G4PrimaryParticle object
+ g4p = new G4PrimaryParticle(IDHEP, PHEP1 * GeV, PHEP2 * GeV, PHEP3 * GeV);
+ g4p->SetMass(PHEP5 * GeV);
+ g4p->SetProperTime(proper_time * ns);
+ _particleMap[p] = g4p;
+ std::set<G4PrimaryParticle*> vec3;
+ for (size_t i = 0; i < p->getDaughters().size(); i++) {
+ if (_visitedParticles.count(p->getDaughters()[i]) == 0) {
+ std::set<G4PrimaryParticle*> vec2 = getRelevantParticles(p->getDaughters()[i]);
+ std::set<G4PrimaryParticle*>::iterator setit;
+ for (setit = vec2.begin(); setit != vec2.end(); setit++) {
+ vec3.insert(*setit);
+ }
+ }
+ }
+ std::set<G4PrimaryParticle*>::iterator setit;
+ for (setit = vec3.begin(); setit != vec3.end(); setit++) {
+ g4p->SetDaughter(*setit);
+ }
+ //std::cout << "*" << std::flush;
+ }
+ //std::cout << g4p->GetPDGcode() << std::flush;
+ relevantParticlesSet.insert(g4p);
+ return relevantParticlesSet;
+ }
+ /* CASE #3:
+ *
+ * If p is not stable AND has zero lifetime: forget about p and retrieve all relevant
+ * decay particles of all daughters of p.
+ *
+ * In this case this step of recursion is just there for collecting all relevant decay
+ * products of daughters and returning them.
+ */
+ else {
+ for (size_t i = 0; i < p->getDaughters().size(); i++) {
+ if (_visitedParticles.count(p->getDaughters()[i]) == 0) {
+ std::set<G4PrimaryParticle*> vec2 = getRelevantParticles(p->getDaughters()[i]);
+ std::set<G4PrimaryParticle*>::iterator setit;
+ for (setit = vec2.begin(); setit != vec2.end(); setit++) {
+ relevantParticlesSet.insert(*setit);
+ }
+ }
+ }
+ return relevantParticlesSet;
+ }
+ }
+ return relevantParticlesSet;
+}
+
+
+} /* namespace slic */
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/StdHepToLcioConvertor.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/StdHepToLcioConvertor.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -1,540 +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"
-
-// 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/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/SteppingAction.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/SteppingAction.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -22,20 +22,20 @@
void SteppingAction::UserSteppingAction(const G4Step* aStep) {
/* TrajectoryManager's stepping action */
- TrajectoryManager::instance()->stepping(aStep);
+ //TrajectoryManager::instance()->stepping(aStep);
/* Check if the track should be killed. */
- checkKillTrack(aStep);
+ //checkKillTrack(aStep);
}
void SteppingAction::checkKillTrack(const G4Step* aStep) {
/* Check if the preStepPoint region has tracking killing enabled. */
- UserRegionInformation* regionInfo = UserRegionInformation::getRegionInformation(aStep->GetPreStepPoint());
- if (regionInfo->getKillTracks()) {
+ //UserRegionInformation* regionInfo = UserRegionInformation::getRegionInformation(aStep->GetPreStepPoint());
+ //if (regionInfo->getKillTracks()) {
/* Kill the track if region track kill flag is set to true. */
- log() << LOG::always << "Killing track: " << aStep->GetTrack()->GetTrackID() << LOG::done;
- aStep->GetTrack()->SetTrackStatus(fStopAndKill);
- }
+ // log() << LOG::always << "Killing track: " << aStep->GetTrack()->GetTrackID() << LOG::done;
+ // aStep->GetTrack()->SetTrackStatus(fStopAndKill);
+ //}
}
} // namespace slic
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/TrackManager.cc (rev 0)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/TrackManager.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -0,0 +1,109 @@
+#include "TrackManager.hh"
+
+// SLIC
+#include "MCParticleManager.hh"
+#include "StdHepGenerator.hh"
+
+// LCIO
+#include "EVENT/LCIO.h"
+#include "EVENT/LCEvent.h"
+
+#include <cassert>
+
+using IMPL::LCCollectionVec;
+using EVENT::LCCollection;
+using EVENT::LCEvent;
+
+namespace slic {
+
+void TrackManager::saveTracks(const G4Event* anEvent, LCEvent* lcEvent) {
+
+ StdHepGenerator::LCIO2Geant4MapIterator it;
+
+ /* Get the map of MCParticles to G4PrimaryParticles. The latter has the track IDs. */
+ StdHepGenerator::LCIO2Geant4Map* particleMap = StdHepGenerator::getParticleMap();
+
+ /* Copy MCParticle information from generated primary particles into the corresponding TrackSummary objects. */
+ for (it = particleMap->begin(); it != particleMap->end(); it++) {
+
+ if (it->second->GetTrackID() < 0) continue;
+
+ for(int j = _trackSummaries->size() - 1; j >= 0; j--){
+ if( _trackSummaries->operator[](j)->getTrackID() == it->second->GetTrackID() ) {
+ _trackSummaries->operator[](j)->setMCParticle(dynamic_cast<MCParticleImpl*>(it->first));
+ break ;
+ }
+ }
+ }
+
+ /* Set parents to be saved also on particles that will be persisted. */
+ for (G4int k = _trackSummaries->size()-1; k >= 0; k--) {
+ if (((*_trackSummaries)[k])->getToBeSaved())
+ ((*_trackSummaries)[k])->setParentToBeSaved();
+ }
+
+ /* Get the MCParticle collection created by event generation. */
+ LCCollectionVec* mcpVec = MCParticleManager::instance()->getMCParticleCollection();
+
+#ifdef SLIC_LOG
+
+#endif
+
+ /* Save TrackSummary objects to LCIO collection. */
+ size_t l;
+ TrackSummary* trackSummary;
+ for (l = 0; l < _trackSummaries->size(); l++) {
+ trackSummary = (*_trackSummaries)[l];
+ if (trackSummary->getToBeSaved()) {
+
+ /* Build the MCParticle in case it hasn't been created yet. */
+ if (trackSummary->getMCParticle() == 0) {
+ trackSummary->buildMCParticle();
+ }
+
+ /* Associate the track ID to the MCParticle for the hit maker. */
+ MCParticleManager::instance()->addMCParticleTrackID(
+ trackSummary->getTrackID(),
+ trackSummary->getMCParticle());
+
+ /*
+ * Only particles created in the simulation need to be added as the primaries
+ * were put into the collection during event generation.
+ */
+ if(trackSummary->getMCParticle()->getGeneratorStatus() == 0 ) {
+ if (trackSummary->getMCParticle() != 0) {
+ mcpVec->push_back(trackSummary->getMCParticle());
+ }
+ }
+ }
+ else {
+ /*
+ * If the track will not be saved then associate the information to its first persisted ancestor.
+ * This will allow association between track IDs from sensitive detector's in LCDD to their appropriate
+ * MCParticle objects in the output collection.
+ */
+ findFirstSavedAncestor(trackSummary);
+ }
+ }
+
+ /* Save the MCParticle collection into the event. */
+ lcEvent->addCollection((LCCollection*)mcpVec, "MCParticle");
+}
+
+TrackSummary* TrackManager::findFirstSavedAncestor(TrackSummary* trackSummary) {
+ TrackSummary* parent = trackSummary->findParent();
+ while (parent != 0) {
+ if (parent->getToBeSaved() == true) {
+ /* Associate the track ID to the first MCParticle up the tree that will be saved. */
+ MCParticleManager::instance()->addMCParticleTrackID(
+ trackSummary->getTrackID(),
+ parent->getMCParticle());
+ break;
+ }
+ parent = parent->findParent();
+ }
+ return parent;
+}
+
+
+}
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/TrackSummary.cc (rev 0)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/TrackSummary.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -0,0 +1,228 @@
+#include "TrackSummary.hh"
+
+// SLIC
+#include "TrackUtil.hh"
+#include "TrackManager.hh"
+
+// Geant4
+#include "G4Track.hh"
+#include "G4DecayProducts.hh"
+#include "G4VProcess.hh"
+
+// LCIO
+#include "EVENT/MCParticle.h"
+
+using EVENT::MCParticle;
+using EVENT::MCParticleVec;
+
+namespace slic {
+
+G4Allocator<TrackSummary> TrackSummaryAllocator;
+
+TrackSummary::TrackSummary(const G4Track* track, G4bool save, G4bool backScattering) :
+ _parent(NULL), _save(save), _backScattering(backScattering), _theMCParticle(NULL),
+ _mcParticleIsUpToDate(false) {
+ _charge = track->GetDefinition()->GetPDGCharge();
+ _mass = track->GetDynamicParticle()->GetMass();
+ _PDG = track->GetDefinition()->GetPDGEncoding();
+ _trackID = track->GetTrackID();
+ _parentTrackID = track->GetParentID();
+}
+
+void TrackSummary::update(const G4Track* track) {
+
+ // TODO: These need to be checked.
+ _momentum = track->GetVertexMomentumDirection();
+ _energy = track->GetTotalEnergy(); // or kinetic E?
+ _vertex = track->GetVertexPosition();
+
+ // update final track information
+ _endPoint = track->GetPosition();
+ _trackLength = track->GetTrackLength();
+
+ // hepEvtStatus == generator status in LCIO
+ _hepEvtStatus = 0;
+ if (_parentTrackID == 0) {
+ const G4DecayProducts* thePreAssignedDecayProducts = track->GetDynamicParticle()->GetPreAssignedDecayProducts();
+ if (thePreAssignedDecayProducts && thePreAssignedDecayProducts->entries() > 0)
+ _hepEvtStatus = 2;
+ else
+ _hepEvtStatus = 1;
+ }
+
+ // Setting the simulator status as for LCIO, by the moment by hand
+ // just adapting the beta LCIO code.
+ // Remarks:
+ // - with Mokka BITVertexIsNotEndpointOfParent is always false
+ // because the endpoint is set always explicitly.
+ // - BITDecayedInTracker and BITDecayedInCalorimeter not
+ // set, waiting for better understand about these bits.
+ //
+ _simstatus = 0;
+
+ // endpoint always set explicitly by Mokka.
+ _simstatus[MCParticle::BITEndpoint] = 1;
+
+ // these folowwing pointers should always be not NULL but
+ // with Geant4 we never know...
+ //
+ const G4Step* theLastStep = track->GetStep();
+ G4StepPoint* theLastPostStepPoint = NULL;
+ if(theLastStep)
+ theLastPostStepPoint = theLastStep->GetPostStepPoint();
+
+ // BITCreatedInSimulation and BITBackscatter has no sense except
+ // if created by simulator.
+ if(_hepEvtStatus == 0)
+ {
+ //fg: BITCreatedInSimulation should only be set if MCParticle has been created
+ // by Mokka/geant4 - now done in BuildMCParticle
+ // _simstatus[ BITCreatedInSimulation ] = 1;
+ if(getBackScattering()) _simstatus[MCParticle::BITBackscatter] = 1;
+ }
+
+ // BITLeftDetector
+ G4bool LeftDetector = false;
+ //
+
+ if(theLastPostStepPoint)
+ //
+ // Thanks to Predrag Krstonosic ([log in to unmask]) we realize that Geant4 never
+ // sets the fWorldBoundary step status even when the particle is leaving the detector
+ // world volume. As a turn around fix while Geant4 does't work correctly we
+ // replaced the test
+ //
+ // LeftDetector = (theLastPostStepPoint->GetStepStatus() == fWorldBoundary))
+ //
+ // by
+ //
+ // LeftDetector = (theLastPostStepPoint->GetStepStatus() == fGeomBoundary))
+ //
+ // in the TrackSummary::update() method. We suppose that normaly a track doens't
+ // finish in a Geometry boundary except if it's leaving the detector world volume,
+ // so it should works nice for the moment.
+ //
+ if((LeftDetector = (theLastPostStepPoint->GetStepStatus() == fGeomBoundary)))
+ _simstatus[MCParticle::BITLeftDetector] = 1;
+
+ // BITStopped
+ if(track->GetKineticEnergy() <= 0.)
+ _simstatus[ MCParticle::BITStopped] = 1;
+
+ G4bool insideTrackerRegion = TrackUtil::getRegionInformation(track)->getStoreSecondaries();
+
+ G4VPhysicalVolume* volume = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->LocateGlobalPointAndSetup(track->GetPosition());
+ G4Region* region = volume->GetLogicalVolume()->GetRegion();
+ UserRegionInformation* regionInfo = dynamic_cast<UserRegionInformation*>(region->GetUserInformation());
+
+ if(insideTrackerRegion)
+ _simstatus[MCParticle::BITDecayedInTracker] = 1;
+
+ // BITDecayedInCalorimeter
+ if(!insideTrackerRegion && !LeftDetector)
+ _simstatus[MCParticle::BITDecayedInCalorimeter] = 1;
+
+ // BITVertexIsNotEndpointOfParent
+ //
+ TrackSummary* myParent = findParent();
+ if(myParent)
+ if(myParent->getEndPoint().isNear(getVertex()))
+ _simstatus[MCParticle::BITVertexIsNotEndpointOfParent] = 1;
+}
+
+void TrackSummary::setParentToBeSaved() {
+ if(!(_theMCParticle && _mcParticleIsUpToDate))
+ buildMCParticle();
+ else
+ // In LCIO mode, if the theMCParticle exists we have already
+ // passed by here.
+ return;
+
+ TrackSummary* myParent = findParent();
+ if (myParent) {
+
+ // std::cout << " **** TrackSummary::SetParentToBeSaved for trackID: " << myParent->GetTrackID() << " this trackid: " << GetTrackID() << std::endl ;
+
+ myParent->_save = true;
+ myParent->setParentToBeSaved();
+
+ MCParticle* theParentMCParticle = myParent->getMCParticle();
+
+ // fg: need to check if particle already has this parent assigned ...
+ // this will be fixed in a future (>1.6) LCIO version
+ const MCParticleVec& parents = _theMCParticle->getParents();
+
+ if(std::find(parents.begin(),
+ parents.end(),
+ theParentMCParticle ) == parents.end() ) {
+
+ // fix B.Vormwald: only add parent if particlke created in simulation
+ if (_theMCParticle->getGeneratorStatus()== 0) {
+ _theMCParticle->addParent( theParentMCParticle );
+ }
+ }
+ }
+}
+
+TrackSummary* TrackSummary::findParent() const {
+
+ TrackSummary* trackSummary = NULL;
+ if (_parentTrackID == 0)
+ return trackSummary;
+
+ TrackSummaryVector* tracks = TrackManager::instance()->getTrackSummaries();
+
+ G4int itr;
+ for (itr = tracks->size() - 1; itr >= 0; itr--) {
+ trackSummary = (*tracks)[itr];
+ if (trackSummary && trackSummary->getTrackID() == _parentTrackID)
+ break;
+ }
+ return trackSummary;
+}
+
+void TrackSummary::buildMCParticle() {
+ // if we have an MCParticle already we just need to update it
+ // otherwise we need to create one ...
+ if( _theMCParticle == 0 ) {
+ _theMCParticle = new MCParticleImpl;
+ _simstatus[MCParticle::BITCreatedInSimulation] = 1; // created in simulation
+ _theMCParticle->setGeneratorStatus(getHepEvtStatus());
+ _theMCParticle->setPDG(getPDG());
+ }
+
+ _theMCParticle->setCharge(getCharge());
+
+ /* Set endpoint. */
+ double endpoint[3];
+ endpoint[0] = getEndPoint()(0);
+ endpoint[1] = getEndPoint()(1);
+ endpoint[2] = getEndPoint()(2);
+ _theMCParticle->setEndpoint(endpoint);
+
+ /* Set mass. */
+ _theMCParticle->setMass(_mass / GeV);
+
+ /* Set simulator status. */
+ _theMCParticle->setSimulatorStatus(getSimulatorStatus());
+
+ /* Set momentum. */
+ float p[3];
+ G4ThreeVector theMomentum = getMomentum();
+ p[0] = theMomentum(0) / GeV;
+ p[1] = theMomentum(1) / GeV;
+ p[2] = theMomentum(2) / GeV;
+ _theMCParticle->setMomentum(p);
+
+ /* Set vertex. */
+ double vertex[3];
+ vertex[0] = getVertex()(0);
+ vertex[1] = getVertex()(1);
+ vertex[2] = getVertex()(2);
+ _theMCParticle->setVertex(vertex);
+
+ /* Set up to date. */
+ _mcParticleIsUpToDate = true;
+}
+
+}
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/TrackingAction.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/TrackingAction.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -1,28 +1,65 @@
// $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/hits/TrackInformation.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() {
-}
+ G4bool isPrimary = (track->GetParentID() == 0);
+ G4ThreeVector trackPosition = track->GetPosition();
+ UserRegionInformation* regionInfo = TrackUtil::getRegionInformation(track);
+ G4bool insideTrackerRegion = regionInfo->getStoreSecondaries();
+ G4bool decayInsideTracker = false;
+ TrackInformation* trackInfo = dynamic_cast<TrackInformation*>(track->GetUserInformation());
-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);
+ /* Track is primary? */
+ if (!isPrimary) {
+
+ /* Inside tracker region? */
+ if (insideTrackerRegion) {
+ const G4VProcess* theCreatorProcess = track->GetCreatorProcess();
+
+ /* Decay in tracker? */
+ if (theCreatorProcess)
+ decayInsideTracker = theCreatorProcess->GetProcessName() == "Decay";
+ }
+ }
+
+ /* Store trajectory? */
+ TrackSummary* trackSummary = 0;
+ if (isPrimary || (insideTrackerRegion && track->GetKineticEnergy() > regionInfo->getThreshold())) {
+
+ /* Turn on trajectory storage in Geant4. */
+ fpTrackingManager->SetStoreTrajectory(true);
+
+ /* Setup the UserTrackInformation if it does not exist. */
+ if (!trackInfo) {
+ TrackUtil::setupUserTrackInformation(track, true);
+ }
+ }
+ /* Turn off trajectory storing. */
+ else {
+ fpTrackingManager->SetStoreTrajectory(false);
+ /* Setup track information, which will not be saved. */
+ TrackUtil::setupUserTrackInformation(track, false);
+ }
}
-void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) {
- // TrajectoryManager's action
- m_trjMgr->postTracking(aTrack);
+void TrackingAction::PostUserTrackingAction(const G4Track* track) {
+ /* Update the TrackSummary with final information about the track. */
+ UserTrackInformation* trackInfo = TrackUtil::getUserTrackInformation(track);
+ if (trackInfo) {
+ trackInfo->getTrackSummary()->update(track);
+ }
}
}
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/TrajectoryManager.cc 2014-01-08 22:29:26 UTC (rev 2863)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/TrajectoryManager.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -147,8 +147,7 @@
nconstTrack->SetUserInformation(trkInfo);
}
-void TrajectoryManager::setTrackingFlagFromRegionInformation(TrackInformation* trk_info,
- UserRegionInformation* reg_info) {
+void TrajectoryManager::setTrackingFlagFromRegionInformation(TrackInformation* trk_info, UserRegionInformation* reg_info) {
if (reg_info->getStoreSecondaries() == true) {
trk_info->setTrackingStatus(TrackInformation::eInTrackingRegion);
} else {
@@ -338,6 +337,12 @@
void TrajectoryManager::stepping(const G4Step* aStep) {
const G4Track* theTrack = aStep->GetTrack();
+ if (inWorldVolume(theTrack)) {
+ std::cout << "stepping in world volume: " << aStep->GetPreStepPoint()->GetPosition() << std::endl;
+ } else {
+ std::cout << "not in world volume: " << aStep->GetPreStepPoint()->GetPosition() << std::endl;
+ }
+
// step has a live track
if (isAlive(aStep)) {
projects/slic/branches/v00-00-01_SLIC-218-dev/src
--- projects/slic/branches/v00-00-01_SLIC-218-dev/src/UserTrackInformation.cc (rev 0)
+++ projects/slic/branches/v00-00-01_SLIC-218-dev/src/UserTrackInformation.cc 2014-01-09 20:37:54 UTC (rev 2864)
@@ -0,0 +1,6 @@
+#include "UserTrackInformation.hh"
+
+namespace slic {
+G4Allocator<UserTrackInformation> TrackInformationAllocator;
+}
+
SVNspam 0.1