slic/src
diff -u -r1.17 -r1.18
--- LcioHitsCollectionBuilder.cc 15 Oct 2009 00:07:36 -0000 1.17
+++ LcioHitsCollectionBuilder.cc 15 Oct 2009 19:14:39 -0000 1.18
@@ -1,4 +1,4 @@
-// $Header: /cvs/lcd/slic/src/LcioHitsCollectionBuilder.cc,v 1.17 2009/10/15 00:07:36 jeremy Exp $
+// $Header: /cvs/lcd/slic/src/LcioHitsCollectionBuilder.cc,v 1.18 2009/10/15 19:14:39 jeremy Exp $
#include "LcioHitsCollectionBuilder.hh"
@@ -344,39 +344,39 @@
// add an MCParticle hit contribution from G4 to LCIO
void LcioHitsCollectionBuilder::addMcpContribs( G4CalorimeterHit* g4CalHit, IMPL::SimCalorimeterHitImpl* simCalHit)
{
- // create empty hit contrib list
+ // Create empty hit contrib list.
McpHitContribList contribs;
- // aggregation of contribs by track ID if CHBIT_PDG is not set
- if ( !m_calCollFlag.bitSet(LCIO::CHBIT_PDG) ) {
-
- // pass ref to contrib list, which will get filled
+ // Use aggregation of contribs by track ID if CHBIT_PDG is not set.
+ if ( !m_calCollFlag.bitSet(LCIO::CHBIT_PDG) )
+ {
+ // Pass a ref to contrib list, which will get filled.
combineMcpHitContribs( g4CalHit->getMcpHitContribList(),
contribs );
}
- // otherwise use the complete list from cal hit
- else {
-
- // set contrib list to cal contribs
+ // Otherwise, use the complete list from the CalHit.
+ else
+ {
contribs = g4CalHit->getMcpHitContribList();
}
- // add contribs to LCIO Mcp
+ // Add contribs to the LCIO MCParticle.
size_t ncontrib = 0;
for ( McpHitContribList::const_iterator iter = contribs.begin();
iter != contribs.end();
- iter++ ) {
+ iter++ )
+ {
- // get this contrib
+ // This contrib.
const McpHitContrib& contrib = ( *iter );
- // fetch associcated MCP pointer by trackID
+ // Get the MCParticle pointer from the track ID.
MCParticleImpl* contribMcp = m_mcpManager->getMaps()->findMcpFromTrackID( contrib.getTrackID() );
- if ( contribMcp != 0 ) {
-
- // add the Mcp contrib to the Lcio cal hit
+ if ( contribMcp != 0 )
+ {
+ // Add the MCParticle contribution to the hit.
simCalHit->addMCParticleContribution(contribMcp,
contrib.getEdep() / GeV,
contrib.getGlobalTime(),
@@ -384,10 +384,18 @@
);
++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;
+ }
+#endif
}
#ifdef SLIC_LOG
- if ( ncontrib == 0 ) {
+ if ( ncontrib == 0 )
+ {
log().error("No hit contribs for sim cal hit.");
}
#endif
slic/src
diff -u -r1.13 -r1.14
--- LcioMcpFactory.cc 15 Oct 2009 00:07:26 -0000 1.13
+++ LcioMcpFactory.cc 15 Oct 2009 19:14:39 -0000 1.14
@@ -1,4 +1,4 @@
-// $Header: /cvs/lcd/slic/src/LcioMcpFactory.cc,v 1.13 2009/10/15 00:07:26 jeremy Exp $
+// $Header: /cvs/lcd/slic/src/LcioMcpFactory.cc,v 1.14 2009/10/15 19:14:39 jeremy Exp $
#include "LcioMcpFactory.hh"
// slic
@@ -44,26 +44,165 @@
// new coll for final Mcps
m_finalColl = static_cast<LCCollectionVec*>(m_manager->getFinalMcpCollection());
- /* Set current G4TrajectoryContainer */
+ // Set current G4TrajectoryContainer
m_currentTrajectoryContainer = event->GetTrajectoryContainer();
- // using StdHep or LCIO e.g. a file-based source?
+ // Check if using StdHep or LCIO e.g. a file-based source?
if ( EventSourceManager::instance()->isFileSource() ) {
-
- // create from the initial coll of Mcp
+ // 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. */
+ // G4 source using trajectories only.
+ else
+ {
+ // Create from trajectory container of current G4Event.
createFinalMcpCollectionFromTrajectoryContainer( m_currentTrajectoryContainer );
}
- /* Fill supplementary collection storing Mcp endpoint energies. */
+ // Fill supplementary collection with MCParticle endpoint energies.
fillMcpEndPointEnergy(m_finalColl);
}
+ IMPL::MCParticleImpl* LcioMcpFactory::createMcpFromTrajectory( Trajectory* trj)
+ {
+#ifdef SLIC_LOG
+ log() << LOG::debug << "createMcpFromTrajectory() - trackId <" << trj->GetTrackID() << ">." << 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 );
+
+ // 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 );
+
+ // time
+ mcp->setTime( trj->GetGlobalTime() );
+
+ // set status codes from trajectory info (with no dep on any associated Mcps)
+ LcioMcpStatusSetter::setMcpStatusCodesFromTrajectory( 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()->addTrackIDToMcpLink( trj->GetTrackID(), mcp );
+
+ return mcp;
+ }
+
+ // FIXME: This method doesn't seem to work correctly!
+ void LcioMcpFactory::createFinalMcpCollectionFromTrajectoryContainer(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 = createMcpFromTrajectory( 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
+ addMcpDaughtersFromTrajectoryContainer( mcp, trj->GetTrackID() );
+ }
+ }
+ }
+ else
+ {
+ G4Exception("G4TrajectoryContainer is null!");
+ }
+
+ // DEBUG
+ LcioMcpManager::instance()->getMaps()->printTrackToMcpMap();
+
+ }
+
+ void LcioMcpFactory::addMcpDaughtersFromTrajectoryContainer( 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()->findMcpFromTrackID( trj->GetTrackID() );
+
+ if ( dauMcp == 0 )
+ {
+ dauMcp = createMcpFromTrajectory( trj );
+ }
+
+ if ( dauMcp == 0 )
+ G4Exception("Failed to create MCParticle.");
+
+ // add this to final coll
+ m_finalColl->addElement( dauMcp );
+
+ // Recursively add daughters of this daughter.
+ addMcpDaughtersFromTrajectoryContainer( dauMcp, trj->GetTrackID() );
+
+ // Set the parent particle in the daughter.
+ dauMcp->addParent( parMcp );
+ }
+ }
+ }
+
+
void LcioMcpFactory::createFinalMcpCollectionFromInitial(EVENT::LCCollection* mcpVecInitial)
{
// coll must exist
@@ -163,61 +302,50 @@
IMPL::MCParticleImpl* LcioMcpFactory::createMcpFromPrimary(G4PrimaryParticle* primary,
EVENT::MCParticle* mcpInit)
{
+#ifdef SLIC_LOG
log() << LOG::done << "createMcpFromPrimary: " << m_manager->getMCParticleIndex( m_manager->getInitialMcpCollection(), mcpInit) << LOG::done;
+#endif
int trkID = primary->GetTrackID();
MCParticleImpl* mcp = 0;
Trajectory* trj = TrajectoryManager::instance()->findTrajectory( trkID );
- /*
- * No trajectory?
- *
- * This is a trap for G4Primaries with unrecognized PDG IDs,
- * aka the dreaded Unknown Particles. In this case, no
- * trajectory was created, so the final Mcp must be composed
- * from initial MCParticle and G4PrimaryParticle information, only.
- *
- */
+ // No trajectory was created.
if ( !trj ) {
#ifdef SLIC_LOG
log() << LOG::debug << "initialAndPrimary" << LOG::done;
#endif
- /* Recursively create Mcps from primary and initial Mcp. */
+ // Recursively create particles from primary and initial MCParticles.
mcp = createMcpFromInitialAndPrimary( primary, mcpInit );
}
- /* Found a trajectory. */
+ // Found a trajectory.
else {
#ifdef SLIC_LOG
log() << LOG::debug << "initialAndTrajectory" << LOG::done;
#endif
- /* Recursively create Mcps from trajectory and initial Mcp. */
+ // Recursively create MCParticles from a trajectory and the initial MCParticle.
mcp = createMcpFromInitialAndTrajectory( trj, mcpInit );
}
- /*
- * In either case, the generator status is set from the initial particle.
- */
+ // Set generator status from 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.
- */
+ // Primary particles do not originate in the simulation.
mcp->setCreatedInSimulation( false );
-
- /*
- * Link the input MCP to an output MCP in order to check later if it was
- * already created.
- */
+
+ // Link the input MCP to an output MCP.
m_manager->getMaps()->addInitialMcpToFinalMcpLink( mcpInit, mcp);
+
+ // Link the track ID to the output particle.
+ m_manager->getMaps()->addTrackIDToMcpLink( trkID, mcp );
+
return mcp;
}
@@ -279,7 +407,7 @@
// found daughter
if ( thisTrj->GetParentID() == trj->GetTrackID() ) {
- log() << LOG::done << "trajDau: " << j << LOG::done;
+// log() << LOG::done << "trajDau: " << j << LOG::done;
// find matching primary Mcp input particle
MCParticle* mcpPrim = m_manager->getMaps()->findPrimaryInitialMcpFromTrajectory( thisTrj );
@@ -346,70 +474,12 @@
// mass
mcp->setMass( primary->GetMass() / GeV );
- /* FIXME: Wrong if G4PrimaryParticle actually has valid charge. (Can it happen here???) */
+ // FIXME: Wrong if G4PrimaryParticle actually has valid charge. (Can it happen?)
mcp->setCharge( LcioMcpManager::m_NAN );
return mcp;
}
- IMPL::MCParticleImpl* LcioMcpFactory::createMcpFromTrajectory( Trajectory* trj)
- {
- log() << LOG::debug << "createMcpFromTrajectory" << LOG::done;
- log() << LOG::debug << "particle: " << trj->GetParticleDefinition()->GetParticleName() << LOG::done;
-
- // 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 );
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "begin pnt: " << beginPos.x() << " " << beginPos.y() << " " << beginPos.z() << LOG::done;
-#endif
-
- // 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 );
-
-#ifdef SLIC_LOG
- log() << LOG::debug << "end pnt: " << endPos.x() << " " << endPos.y() << " " << endPos.z() << LOG::done;
-#endif
-
- // time
- mcp->setTime( trj->GetGlobalTime() );
-
- // add to map of trkID -> mcp
- m_manager->getMaps()->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,
@@ -447,75 +517,7 @@
primDau = primDau->GetNext();
}
}
-
- void LcioMcpFactory::createFinalMcpCollectionFromTrajectoryContainer(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 ) {
-
- // 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
- 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()->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)
{