Author: [log in to unmask] Date: Tue Oct 13 14:16:49 2015 New Revision: 3684 Log: [SLIC-252] Make sure particles with gen status 0 are correctly handled; add additional sanity checks; minor change to debug output. Modified: projects/slic/trunk/src/MCParticleManager.cc Modified: projects/slic/trunk/src/MCParticleManager.cc ============================================================================= --- projects/slic/trunk/src/MCParticleManager.cc (original) +++ projects/slic/trunk/src/MCParticleManager.cc Tue Oct 13 14:16:49 2015 @@ -13,7 +13,7 @@ void MCParticleManager::generateEvent(LCCollectionVec* particles, G4Event* event) { #if SLIC_LOG - log() << LOG::head << LOG::debug << "generating event " << event->GetEventID() << " from " << particles->size() << " MCParticles" << LOG::done; + log() << LOG::debug << LOG::head << "generating event " << event->GetEventID() << " from " << particles->size() << " particles" << LOG::endl << LOG::done; #endif // Apply Z smearing to input particles. @@ -38,25 +38,23 @@ // Debug print information about the particle. #ifdef SLIC_LOG - log() << LOG::debug << "processing MCParticle " << i << LOG::done; - log() << LOG::debug - << "PDG: " << mcp->getPDG() << LOG::endl - << "gen status: " << mcp->getGeneratorStatus() << LOG::endl - << "sim status: " << mcp->getSimulatorStatus() << LOG::endl - << "energy: " << mcp->getEnergy() << LOG::endl - << "time: " << mcp->getTime() << LOG::endl - << "daughters: " << mcp->getDaughters().size() << LOG::endl - << "vertex: ( " << mcp->getVertex()[0] << " " << mcp->getVertex()[1] << " " << mcp->getVertex()[2] << " )" << LOG::endl - << "p: ( " << mcp->getMomentum()[0] << ", " << mcp->getMomentum()[1] << ", " << mcp->getMomentum()[2] << " )" << LOG::endl - << "charge: " << mcp->getCharge() << LOG::endl - << "mass: " << mcp->getMass() << LOG::endl - << "daughters: " << mcp->getDaughters().size() << LOG::endl - << "parents: " << mcp->getParents().size() + log() << LOG::debug << "processing particle " << i << LOG::endl + << " PDG: " << mcp->getPDG() << LOG::endl + << " gen status: " << mcp->getGeneratorStatus() << LOG::endl + << " sim status: " << mcp->getSimulatorStatus() << LOG::endl + << " energy: " << mcp->getEnergy() << LOG::endl + << " time: " << mcp->getTime() << LOG::endl + << " daughters: " << mcp->getDaughters().size() << LOG::endl + << " vertex: ( " << mcp->getVertex()[0] << " " << mcp->getVertex()[1] << " " << mcp->getVertex()[2] << " )" << LOG::endl + << " p: ( " << mcp->getMomentum()[0] << ", " << mcp->getMomentum()[1] << ", " << mcp->getMomentum()[2] << " )" << LOG::endl + << " charge: " << mcp->getCharge() << LOG::endl + << " mass: " << mcp->getMass() << LOG::endl + << " daughters: " << mcp->getDaughters().size() << LOG::endl + << " parents: " << mcp->getParents().size() << LOG::endl << LOG::done; - log() << LOG::endl; -#endif - - // Only make vertices for particles without parents. +#endif + + // Only process particles without parents; daughters will be handled recursively. if (mcp->getParents().size() == 0) { #ifdef SLIC_LOG @@ -64,41 +62,65 @@ #endif // Create a new primary vertex. - G4PrimaryVertex* vertex = 0; - - if (mcp->getGeneratorStatus() == 1 || mcp->getGeneratorStatus() == 2) { - vertex = createVertex(mcp); + G4PrimaryVertex* vertex = createVertex(mcp); + + if (vertex == 0) { +#ifdef SLIC_LOG + log() << LOG::debug << "no vertex created for particle " << i << LOG::done; +#endif } #ifdef SLIC_LOG log() << LOG::debug << "creating primary for MCParticle " << i << LOG::done; #endif - // Create a primary particle and recursively create its daughters. + // Create a primary particle and recursively create primaries for its daughters. const std::set<G4PrimaryParticle*> primaries = createPrimary(0, mcp); - // Loop over returned set of primary particles. - for (std::set<G4PrimaryParticle*>::const_iterator it = primaries.begin(); it != primaries.end(); it++) { - - // Is this a top-level particle with no daughters? - if (!MCParticleManager::isDaughter((*it), primaries)) { - - // Create a new vertex if necessary. - if (vertex == 0) { - // Create vertex by looking up MCParticle from the primary. - vertex = createVertex(_mcpMap[(*it)]); + // Were any primaries created? + if (primaries.size() > 0) { + + // Loop over returned set of primary particles. + for (std::set<G4PrimaryParticle*>::const_iterator it = primaries.begin(); it != primaries.end(); it++) { + + // Is this a top-level particle with no daughters? + if (!MCParticleManager::isDaughter((*it), primaries)) { + + // Create a new vertex if necessary. + if (vertex == 0) { +#ifdef SLIC_LOG + log() << "creating vertex from primary " << *it << LOG::done; +#endif + // Create vertex by looking up MCParticle from the primary. + vertex = createVertex(_mcpMap[(*it)]); + } + + // Add the primary particle to the vertex. +#ifdef SLIC_LOG + log() << "adding primary daughter particle " << *it << " to vertex" << LOG::done; +#endif + vertex->SetPrimary(*it); } - - // Add the primary particle to the vertex. - vertex->SetPrimary(*it); + } + + // Add the vertex to the event. + event->AddPrimaryVertex(vertex); + + } else { +#ifdef SLIC_LOG + log() << LOG::warning << "skipping particle " << i << " with no primaries" << LOG::done; +#endif + // No primaries created so this record and its daughters are skipped! + if (vertex != 0) { +#ifdef SLIC_LOG + log() << LOG::debug << "deleting empty vertex for particle " << i << LOG::done; +#endif + delete vertex; } } - - // Add the vertex to the event. - event->AddPrimaryVertex(vertex); - } -#ifdef SLIC_LOG - log() << "-------------------------------------------------" << LOG::endl; + } +#ifdef SLIC_LOG + log() << LOG::debug << LOG::endl << "-------------------------------------------------" << LOG::endl << LOG::done; #endif } @@ -141,7 +163,7 @@ // Map primary to MCParticle. #ifdef SLIC_LOG - log() << LOG::debug << "mapping MCParticle " << mcp << " to primary " << primary << LOG::done; + log() << LOG::debug << "mapping particle " << mcp << " to primary " << primary << LOG::done; #endif _primaryMap[mcp] = primary; @@ -160,7 +182,7 @@ for (int i = 0; i < mcp->getDaughters().size(); i++) { #ifdef SLIC_LOG - log() << LOG::debug << "creating primary for dau particle " << i << LOG::done; + log() << LOG::debug << "creating primary for daughter " << i << LOG::done; #endif std::set<G4PrimaryParticle*> daughters = createPrimary(primary, mcp->getDaughters()[i]); #ifdef SLIC_LOG @@ -177,18 +199,23 @@ G4PrimaryVertex* MCParticleManager::createVertex(MCParticle* mcp) { - // Get the particle's origin. - const double* particleVertex = mcp->getVertex(); - G4ThreeVector particlePosition(particleVertex[0], particleVertex[1], particleVertex[2]); - - // Get the particle's time. - G4double particleTime = mcp->getTime(); - - // Create a new vertex for this particle. -#ifdef SLIC_LOG - log() << LOG::debug << "creating vertex at position " << particlePosition << " with time " << particleTime << LOG::done; -#endif - G4PrimaryVertex* vertex = new G4PrimaryVertex(particlePosition, particleTime); + G4PrimaryVertex* vertex = 0; + + if (mcp->getGeneratorStatus() == 1 || mcp->getGeneratorStatus() == 2) { + + // Get the particle's origin. + const double* particleVertex = mcp->getVertex(); + G4ThreeVector particlePosition(particleVertex[0], particleVertex[1], particleVertex[2]); + + // Get the particle's time. + G4double particleTime = mcp->getTime(); + + // Create a new vertex for this particle. +#ifdef SLIC_LOG + log() << LOG::debug << "creating vertex at position " << particlePosition << " with time " << particleTime << LOG::done; +#endif + vertex = new G4PrimaryVertex(particlePosition, particleTime); + } return vertex; } ######################################################################## Use REPLY-ALL to reply to list To unsubscribe from the LCDET-SVN list, click the following link: https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCDET-SVN&A=1