slic/include
diff -N LcioMcpFactory.hh
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ LcioMcpFactory.hh 14 Sep 2005 23:09:48 -0000 1.1
@@ -0,0 +1,88 @@
+// $Header: /cvs/lcd/slic/include/LcioMcpFactory.hh,v 1.1 2005/09/14 23:09:48 jeremy Exp $
+#ifndef LcioMcpFactory_hh
+#define LcioMcpFactory_hh 1
+
+// LCIO
+#include "EVENT/LCCollection.h"
+#include "EVENT/LCEvent.h"
+#include "EVENT/MCParticle.h"
+#include "IMPL/MCParticleImpl.h"
+
+// slic
+#include "Verbose.hh"
+
+class G4PrimaryParticle;
+class G4TrajectoryContainer;
+
+namespace slic
+{
+ class LcioMcpManager;
+ class Trajectory;
+
+ /**
+ * @class LcioMcpFactory
+ * @brief Creates LCIO MCParticles from saved Geant4/SLIC runtime information.
+ */
+ class LcioMcpFactory : public Verbose
+ {
+
+ public:
+ LcioMcpFactory(LcioMcpManager*);
+ virtual ~LcioMcpFactory();
+
+ public:
+
+ // create new, blank Mcp coll
+ EVENT::LCCollection* createEmptyMcpCollection(const std::string& collName, bool errorOnExist = true);
+
+ // chooses to create from initial or traj cont
+ void createFinalMcpCollection();
+
+ // use StdHep-generated mcpVec to create final mcpVec for an event
+ void createFinalMcpCollectionFromInitial(EVENT::LCCollection* mcpVecInitial);
+
+ // These two functions create the Mcp coll from trajectories only in case of G4 GPS or gun
+ void createFinalMcpCollectionFromTrajectoryContainer(G4TrajectoryContainer* trjCont);
+
+ /* Create Mcp daughters of a trajectory with the given track ID. */
+ void addMcpDaughtersFromTrajectoryContainer(IMPL::MCParticleImpl* parMcp, int parTrkID);
+
+ // create an Mcp from a trajectory only
+ IMPL::MCParticleImpl* createMcpFromTrajectory( Trajectory* );
+
+ // create Mcps recursively given an initial Mcp
+ IMPL::MCParticleImpl* createMcpFromInitialRecurse(EVENT::MCParticle* mcp);
+
+ // create an Mcp from another, associated initial Mcp only
+ IMPL::MCParticleImpl* createMcpFromInitialOnly(EVENT::MCParticle* mcpInit);
+
+ // create an Mcp from a G4PrimaryParticle
+ IMPL::MCParticleImpl* createMcpFromPrimary(G4PrimaryParticle*, EVENT::MCParticle* );
+
+ // create an Mcp from primary and input Mcp
+ IMPL::MCParticleImpl* createMcpFromInitialAndPrimary(G4PrimaryParticle* primary, EVENT::MCParticle* mcpInit);
+
+ // create an Mcp from a trajectory and an input Mcp
+ IMPL::MCParticleImpl* createMcpFromInitialAndTrajectory( Trajectory*, EVENT::MCParticle* mcpInit );
+
+ // create an Mcp from a G4PrimaryParticle, shallow copy only
+ IMPL::MCParticleImpl* createMcpFromPrimaryShallowCopy(G4PrimaryParticle*);
+
+ // create new Mcp and shallow copy from input Mcp
+ IMPL::MCParticleImpl* createMcpShallowCopy(EVENT::MCParticle* mcp);
+
+ // create and add MCP daughters based on G4PrimaryParticle and associated MCP from StdHep
+ void createDaughtersFromPrimary(G4PrimaryParticle* primary,
+ EVENT::MCParticle* mcpInit,
+ IMPL::MCParticleImpl* mcpPar);
+
+ // add daughters to a MCP based on associated intial MCP from StdHep
+ void addMcpDaughtersFromInitial( IMPL::MCParticleImpl* mcpNew, EVENT::MCParticle* mcpInit );
+
+ private:
+ LcioMcpManager* m_manager;
+ EVENT::LCCollection* m_finalColl;
+ };
+}
+
+#endif
slic/src
diff -N LcioMcpFactory.cc
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ LcioMcpFactory.cc 14 Sep 2005 23:09:49 -0000 1.1
@@ -0,0 +1,472 @@
+// $Header: /cvs/lcd/slic/src/LcioMcpFactory.cc,v 1.1 2005/09/14 23:09:49 jeremy Exp $
+#include "LcioMcpFactory.hh"
+
+// slic
+#include "LcioMcpManager.hh"
+#include "LcioMcpUtil.hh"
+#include "Trajectory.hh"
+#include "EventSourceManager.hh"
+#include "TrajectoryManager.hh"
+#include "LcioMcpStatusSetter.hh"
+
+// geant4
+#include "globals.hh"
+#include "G4TrajectoryContainer.hh"
+#include "G4EventManager.hh"
+
+// lcio
+#include "IMPL/LCCollectionVec.h"
+
+using EVENT::LCIO;
+using EVENT::MCParticle;
+using IMPL::MCParticleImpl;
+using IMPL::LCCollectionVec;
+
+namespace slic
+{
+
+ LcioMcpFactory::LcioMcpFactory(LcioMcpManager* manager)
+ {
+ m_manager = manager;
+ }
+
+ LcioMcpFactory::~LcioMcpFactory()
+ {}
+
+ void LcioMcpFactory::createFinalMcpCollection()
+ {
+ // new coll for final Mcps
+ m_finalColl = m_manager->getFinalMcpCollection();
+
+ // using StdHep or LCIO e.g. a file-based source?
+ if ( EventSourceManager::instance()->isFileSource() ) {
+
+ // create from the initial coll of Mcp
+ createFinalMcpCollectionFromInitial( m_manager->getInitialMcpCollection() );
+ }
+ // G4 source using traj only
+ else {
+
+ /* Create from trj cont of current G4Event, only. */
+ createFinalMcpCollectionFromTrajectoryContainer(G4EventManager::GetEventManager()
+ ->GetNonconstCurrentEvent()
+ ->GetTrajectoryContainer() );
+ }
+ }
+
+ void LcioMcpFactory::createFinalMcpCollectionFromInitial(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++ ) {
+
+ // get Mcp
+ MCParticleImpl* mcp = static_cast<MCParticleImpl*> ( mcpVecInitial->getElementAt( i ) );
+
+ // create only from Mcp primaries, e.g. no parents
+ if ( LcioMcpUtil::isPrimary( mcp ) ) {
+ createMcpFromInitialRecurse( mcp );
+ }
+ }
+ }
+ else {
+ G4Exception("Initial McpVec has no members.");
+ }
+ }
+ else {
+ G4Exception("Initial McpVec ptr is null.");
+ }
+ }
+
+ IMPL::MCParticleImpl* LcioMcpFactory::createMcpFromInitialRecurse(EVENT::MCParticle* mcpInit)
+ {
+ MCParticleImpl* mcp = 0;
+
+ // find associated primary
+ G4PrimaryParticle* g4primary = m_manager->findPrimaryFromMcp( mcpInit );
+
+ // no associated primary with this initial mcp?
+ if( !g4primary ) {
+
+ // create from initial Mcp only
+ mcp = createMcpFromInitialOnly( mcpInit );
+ }
+ // create from primary
+ else {
+
+ // create it from primary
+ mcp = createMcpFromPrimary( g4primary, mcpInit );
+ }
+
+ // return the output MCP
+ return mcp;
+ }
+
+ IMPL::MCParticleImpl* LcioMcpFactory::createMcpFromInitialOnly(EVENT::MCParticle* mcpInit)
+ {
+ // copy shallow the initial
+ MCParticleImpl* mcp = createMcpShallowCopy( mcpInit );
+
+ // should add to Mcp map here?
+ m_manager->addInitialMcpToFinalMcpLink( mcpInit, mcp);
+
+ // given initial, create and add daughters
+ addMcpDaughtersFromInitial( mcp, mcpInit );
+
+ // add to final coll
+ m_finalColl->addElement( mcp );
+
+ return mcp;
+ }
+
+ IMPL::MCParticleImpl* LcioMcpFactory::createMcpFromPrimary(G4PrimaryParticle* primary,
+ EVENT::MCParticle* mcpInit)
+ {
+ int trkID = primary->GetTrackID();
+ MCParticleImpl* mcp = 0;
+ Trajectory* trj = TrajectoryManager::instance()->findTrajectory( trkID );
+
+ /*
+ * No trajectory?
+ *
+ * This is a trap for G4Primaries with unrecognized PDGid,
+ * aka the dreaded "unknown particle". In this case, no
+ * trajectory was created, so the final Mcp must be composed
+ * from trajectory and initial Mcp information.
+ *
+ */
+ if ( !trj ) {
+
+ /* Recursively create Mcps from primary and initial Mcp. */
+ mcp = createMcpFromInitialAndPrimary( primary, mcpInit );
+
+ }
+ /* Found a trajectory. */
+ else {
+
+ /* Recursively create Mcps from trajectory and initial Mcp. */
+ mcp = createMcpFromInitialAndTrajectory( trj, mcpInit );
+
+ }
+
+ /*
+ * In either case, the generator status is set from the initial particle.
+ */
+ mcp->setGeneratorStatus( mcpInit->getGeneratorStatus() );
+
+ /*
+ * The particle did not originate in the simulator, because it
+ * has a matching primary, which means that it had a corresponding
+ * record from the event generator.
+ */
+ mcp->setCreatedInSimulation( false );
+
+ /*
+ * Link the input MCP to an output MCP in order to check later if it was
+ * already created.
+ */
+ m_manager->addInitialMcpToFinalMcpLink( mcpInit, mcp);
+
+ return mcp;
+ }
+
+ IMPL::MCParticleImpl* LcioMcpFactory::createMcpFromInitialAndPrimary(G4PrimaryParticle* primary,
+ EVENT::MCParticle* mcpInit)
+ {
+ /*
+ * Create new Mcp from primary particle,
+ * shallow copy only (no recursion).
+ */
+ MCParticleImpl* mcp = createMcpFromPrimaryShallowCopy( 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] };
+ mcp->setVertex( vtx );
+
+ /* setup primary's daughters */
+ createDaughtersFromPrimary( primary, /* this primary */
+ mcpInit, /* associated initial Mcp particle */
+ mcp /* parent Mcp particle */
+ );
+
+#ifdef SLIC_VERBOSE
+ if ( mcp->getDaughters().size() == 0 ) {
+ if ( verbose() ) {
+ std::cout << "No Mcp daughters added." << std::endl;
+ }
+ }
+#endif
+
+ // done with particle; add to coll
+ m_finalColl->addElement( mcp );
+
+ return mcp;
+ }
+
+ IMPL::MCParticleImpl* LcioMcpFactory::createMcpFromInitialAndTrajectory(Trajectory* trj,
+ EVENT::MCParticle* mcpInit)
+ {
+ // shallow copy create with basic info and no daughters
+ IMPL::MCParticleImpl* mcp = createMcpFromTrajectory( trj );
+
+ // set generator status and created in sim flag
+ LcioMcpStatusSetter::setGeneratorStatus( mcpInit, mcp );
+
+ // loop over all trajectories to find trajectory daughter
+ G4TrajectoryContainer* trjCont = G4EventManager::GetEventManager()->GetNonconstCurrentEvent()->GetTrajectoryContainer();
+ int numTrj = trjCont->entries();
+ for( int j=0; j < numTrj; j++ ) {
+
+ Trajectory* thisTrj = static_cast<Trajectory*> ((*trjCont)[j]);
+
+ // found daughter
+ if ( thisTrj->GetParentID() == trj->GetTrackID() ) {
+
+ // find matching primary Mcp input particle
+ MCParticle* mcpPrim = m_manager->findPrimaryInitialMcpFromTrajectory( thisTrj );
+
+ // create MCP output particle from this input daughter
+ MCParticleImpl* mcpDauNew = createMcpFromInitialAndTrajectory( thisTrj, mcpPrim );
+
+ // set parent link of new daughter
+ mcpDauNew->addParent( mcp );
+ }
+ }
+
+ // add to collection
+ m_finalColl->addElement( mcp );
+
+ return mcp;
+ }
+
+ IMPL::MCParticleImpl* LcioMcpFactory::createMcpShallowCopy(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() );
+
+ return mcpNew;
+ }
+
+ IMPL::MCParticleImpl* LcioMcpFactory::createMcpFromPrimaryShallowCopy(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 );
+
+ // charge
+ mcp->setCharge( primary->GetCharge() );
+
+ return mcp;
+ }
+
+ IMPL::MCParticleImpl* LcioMcpFactory::createMcpFromTrajectory( Trajectory* trj)
+ {
+ // 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 );
+
+ // mass
+ mcp->setMass( trj->GetParticleDefinition()->GetPDGMass() / GeV );
+
+ // 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 );
+
+ // add to map of trkID -> mcp
+ m_manager->addTrackIDToMcpLink( trj->GetTrackID(), mcp );
+
+ // set status codes from trajectory info (with no dep on any associated Mcps)
+ LcioMcpStatusSetter::setMcpStatusCodesFromTrajectory( trj, mcp );
+
+ return mcp;
+ }
+
+ // create LCIO MCParticles from daughters, setting proper linkage
+ void LcioMcpFactory::createDaughtersFromPrimary(G4PrimaryParticle* primary,
+ EVENT::MCParticle* mcpInit,
+ IMPL::MCParticleImpl* mcpPar)
+ {
+
+#ifdef SLIC_VERBOSE
+ if ( verbose() > 1 ) {
+ std::cout << "LcioMcpFactory::createDaughtersFromPrimary()" << std::endl;
+ std::cout << "primary <" << primary << ">" << std::endl;
+ std::cout << "trkID <" << primary->GetTrackID() << ">" << std::endl;
+ std::cout << "mcpInit <" << mcpInit << ">" << std::endl;
+ std::cout << "mcpPar <" << mcpPar << ">" << std::endl;
+ }
+#endif
+
+ // loop over primary daughters
+ G4PrimaryParticle* primDau = primary->GetDaughter();
+ while ( primDau ) {
+
+ // find input Mcp for the primary dau
+ MCParticle* mcpDau = m_manager->findDaughterMcpFromPrimary( mcpInit,
+ primDau );
+
+ // must find it or die
+ if ( 0 == mcpDau ) {
+ G4Exception( "Mcp daughter was not found." );
+ }
+
+ // create dau Mcp
+ MCParticleImpl* mcpDauNew = createMcpFromPrimary( primDau, mcpDau );
+
+ // make mcpPar the parent
+ mcpDauNew->addParent( mcpPar );
+
+ // get next primary dau
+ primDau = primDau->GetNext();
+ }
+ }
+
+ void LcioMcpFactory::createFinalMcpCollectionFromTrajectoryContainer(G4TrajectoryContainer* trjCont)
+ {
+ if ( trjCont ) {
+
+ // loop over trj cont
+ int n_trj = trjCont->entries();
+ for ( int i = 0;
+ i < n_trj;
+ i++ ) {
+
+ Trajectory* trj = static_cast<Trajectory*> ((*trjCont)[i]);
+
+ // only add primary trajectories; daughters are picked up within
+ if ( trj->GetParentID() == 0 ) {
+
+ // create mcp output particle from trj
+ MCParticleImpl* mcp = createMcpFromTrajectory( trj );
+
+ // add this to coll
+ m_finalColl->addElement( mcp );
+
+ // add daughters
+ addMcpDaughtersFromTrajectoryContainer( mcp, trj->GetTrackID() );
+ }
+ }
+ }
+ else {
+ G4Exception("G4TrajectoryContainer is null!");
+ }
+ }
+
+ void LcioMcpFactory::addMcpDaughtersFromTrajectoryContainer( MCParticleImpl* parMcp,
+ int parTrkID)
+ {
+ // loop over trj cont
+ G4TrajectoryContainer* trjCont = G4EventManager::GetEventManager()->GetNonconstCurrentEvent()->GetTrajectoryContainer();
+ int n_trj = trjCont->entries();
+ for ( int i = 0;
+ i < n_trj;
+ i++ ) {
+
+ Trajectory* trj = static_cast<Trajectory*> ((*trjCont)[i]);
+
+ if ( trj->GetParentID() == parTrkID ) {
+
+ // lkp dau Mcp by track ID
+ MCParticleImpl* dauMcp = m_manager->findMcpFromTrackID( trj->GetTrackID() );
+
+ if ( dauMcp == 0 ) {
+
+ // create from trj if null
+ dauMcp = createMcpFromTrajectory( trj );
+ }
+
+ // must be created by now
+ assert( dauMcp );
+
+ // add this to final coll
+ m_finalColl->addElement( dauMcp );
+
+ // add daughters
+ addMcpDaughtersFromTrajectoryContainer( dauMcp, trj->GetTrackID() );
+
+ // set parent in daughter
+ dauMcp->addParent( parMcp );
+ }
+ }
+ }
+
+ void LcioMcpFactory::addMcpDaughtersFromInitial(IMPL::MCParticleImpl* mcpNew,
+ EVENT::MCParticle* mcpInit)
+ {
+ 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->findFinalParticleFromInitial( mcpChildInit );
+
+ // if not added already
+ if ( mcpChildFinal == 0 ) {
+
+ // create daughter
+ mcpChildFinal = createMcpFromInitialRecurse( mcpChildInit );
+ }
+
+ // add parent regardless of above; maybe new MCP or already created
+ mcpChildFinal->addParent( mcpNew );
+ }
+ }
+}
