Skip to content
Snippets Groups Projects
Commit 16a4b1af authored by Richard Hawkings's avatar Richard Hawkings Committed by Graeme Stewart
Browse files

add merge of pileup event collection in GenEventStackFiller (ISF_HepMC_Tools-00-03-08)

	* src/GenEventStackFiller.cxx, .h
	- add ability to merge in separate McEventCollection for pileup
	  triggered by new parameter PileupMCEventCollection, default is ""
	  which keeps behaviour as before
	* tag as ISF_HepMC_Tools-00-03-08
parent 900fa1e8
No related branches found
No related tags found
No related merge requests found
...@@ -48,6 +48,8 @@ ISF::GenEventStackFiller::GenEventStackFiller(const std::string& t, const std::s ...@@ -48,6 +48,8 @@ ISF::GenEventStackFiller::GenEventStackFiller(const std::string& t, const std::s
m_particleDataTable(0), m_particleDataTable(0),
m_inputMcEventCollection("GEN_EVENT"), m_inputMcEventCollection("GEN_EVENT"),
m_outputMcEventCollection("TruthEvent"), m_outputMcEventCollection("TruthEvent"),
m_pileupMcEventCollection(""),
m_outputPileupMcEventCollection(""),
m_recordOnlySimulated(false), m_recordOnlySimulated(false),
m_useGeneratedParticleMass(false), m_useGeneratedParticleMass(false),
m_genEventManipulators(), m_genEventManipulators(),
...@@ -73,6 +75,11 @@ ISF::GenEventStackFiller::GenEventStackFiller(const std::string& t, const std::s ...@@ -73,6 +75,11 @@ ISF::GenEventStackFiller::GenEventStackFiller(const std::string& t, const std::s
declareProperty("PurgeOutputCollectionToSimulatedParticlesOnly", declareProperty("PurgeOutputCollectionToSimulatedParticlesOnly",
m_recordOnlySimulated, m_recordOnlySimulated,
"Record only particles that were taken as Input for the Simulation."); "Record only particles that were taken as Input for the Simulation.");
declareProperty("PileupMcEventCollection",m_pileupMcEventCollection,
"StoreGate collection name of pileup generator McEventCollection");
declareProperty("OutputPileupMcEventCollection",
m_outputPileupMcEventCollection,
"StoreGate collection name of pileup truth McEventCollection");
// particle mass from particle data table? // particle mass from particle data table?
declareProperty("UseGeneratedParticleMass", declareProperty("UseGeneratedParticleMass",
m_useGeneratedParticleMass, m_useGeneratedParticleMass,
...@@ -164,6 +171,9 @@ StatusCode ISF::GenEventStackFiller::fillStack(ISF::ISFParticleContainer& partic ...@@ -164,6 +171,9 @@ StatusCode ISF::GenEventStackFiller::fillStack(ISF::ISFParticleContainer& partic
} else { } else {
// 2. try to retrieve and process regular mceventcollection, as usual // 2. try to retrieve and process regular mceventcollection, as usual
sc = processSingleColl(particleColl); sc = processSingleColl(particleColl);
// added by RJH, process additional inline pileup collections if requested
if (sc.isSuccess() && m_pileupMcEventCollection!="")
sc = processPileupColl(particleColl);
} }
if (sc.isFailure()) { if (sc.isFailure()) {
...@@ -180,7 +190,7 @@ StatusCode ISF::GenEventStackFiller::fillStack(ISF::ISFParticleContainer& partic ...@@ -180,7 +190,7 @@ StatusCode ISF::GenEventStackFiller::fillStack(ISF::ISFParticleContainer& partic
} }
/** process merged (hard-scatter+pileup) mcevent collections */ /** process single (hard-scatter-only) mcevent collection */
StatusCode StatusCode
ISF::GenEventStackFiller::processSingleColl(ISF::ISFParticleContainer& particleColl) const ISF::GenEventStackFiller::processSingleColl(ISF::ISFParticleContainer& particleColl) const
{ {
...@@ -206,13 +216,13 @@ ISF::GenEventStackFiller::processSingleColl(ISF::ISFParticleContainer& particleC ...@@ -206,13 +216,13 @@ ISF::GenEventStackFiller::processSingleColl(ISF::ISFParticleContainer& particleC
return StatusCode::FAILURE; return StatusCode::FAILURE;
} }
ATH_MSG_VERBOSE( "Created particle collection with length " << particleColl.size() ); ATH_MSG_DEBUG( "Created particle collection with length " << particleColl.size() );
return StatusCode::SUCCESS; return StatusCode::SUCCESS;
} }
/** process single (hard-scatter-only) mcevent collection */ /** process merged (hard-scatter+pileup) mcevent collections */
StatusCode StatusCode
ISF::GenEventStackFiller::processMergedColls(ISF::ISFParticleContainer& particleColl) const ISF::GenEventStackFiller::processMergedColls(ISF::ISFParticleContainer& particleColl) const
{ {
...@@ -264,6 +274,37 @@ ISF::GenEventStackFiller::processMergedColls(ISF::ISFParticleContainer& particle ...@@ -264,6 +274,37 @@ ISF::GenEventStackFiller::processMergedColls(ISF::ISFParticleContainer& particle
return StatusCode::SUCCESS; return StatusCode::SUCCESS;
} }
// addd by RJH process pileup mcevent collection
StatusCode
ISF::GenEventStackFiller::processPileupColl(ISF::ISFParticleContainer&
particleColl) const
{
// Retrieve pileup collecton from StoreGate
const McEventCollection* inputCollection=0;
if (evtStore()->retrieve(inputCollection,m_pileupMcEventCollection).isFailure()) {
ATH_MSG_ERROR("Coould not retrieve " << m_pileupMcEventCollection <<
"from StoreGateSvc");
return StatusCode::FAILURE;
}
// create copy
McEventCollection *mcCollection=new McEventCollection(*inputCollection);
StatusCode status=mcEventCollLooper(particleColl,mcCollection,1);
if (status!=StatusCode::SUCCESS) return status;
// record this additional collection to Storegate
bool allowMods(true);
if (evtStore()->record(mcCollection, m_outputPileupMcEventCollection,
allowMods).isFailure()) {
ATH_MSG_ERROR( "Could not record " << m_outputPileupMcEventCollection <<
" to StoreGateSvc. Abort.");
return StatusCode::FAILURE;
}
ATH_MSG_DEBUG( "After adding pileup, particle collection has length " <<
particleColl.size() );
return StatusCode::SUCCESS;
}
/** common code to loop over mcevent collection */ /** common code to loop over mcevent collection */
StatusCode StatusCode
...@@ -272,6 +313,10 @@ ISF::GenEventStackFiller::mcEventCollLooper(ISF::ISFParticleContainer& particleC ...@@ -272,6 +313,10 @@ ISF::GenEventStackFiller::mcEventCollLooper(ISF::ISFParticleContainer& particleC
// ensure unique barcode for merged pileup collisions // ensure unique barcode for merged pileup collisions
int bcPileupOffset = ( nPileupCounter>0 ? m_largestBc+1 : 0 ); int bcPileupOffset = ( nPileupCounter>0 ? m_largestBc+1 : 0 );
ATH_MSG_DEBUG("Starting mcEVentCollLooper for pileup " << nPileupCounter
<< " bcPileupOffset " << bcPileupOffset << " with "
<< particleColl.size() << " ISF particles");
McEventCollection::iterator eventIt = mcCollection->begin(); McEventCollection::iterator eventIt = mcCollection->begin();
McEventCollection::iterator eventItEnd = mcCollection->end(); McEventCollection::iterator eventItEnd = mcCollection->end();
...@@ -284,10 +329,18 @@ ISF::GenEventStackFiller::mcEventCollLooper(ISF::ISFParticleContainer& particleC ...@@ -284,10 +329,18 @@ ISF::GenEventStackFiller::mcEventCollLooper(ISF::ISFParticleContainer& particleC
std::vector<HepMC::GenParticle*> goodparticles; std::vector<HepMC::GenParticle*> goodparticles;
// loop over the event in the mc collection // loop over the event in the mc collection
for ( ; eventIt != eventItEnd; ++eventIt) { int npile=0;
for ( ; eventIt != eventItEnd; ++eventIt,++npile) {
// skip empty events // skip empty events
if ( *eventIt == 0) continue; if ( *eventIt == 0) continue;
// RJH - for pileup, calculate bcid from signal-process-id
int bcid=0;
if (nPileupCounter>0) {
bcid=((*eventIt)->signal_process_id())/10000;
ATH_MSG_DEBUG("Pileup index " << npile << " mapped to BCID " << bcid);
}
ATH_MSG_VERBOSE(" signal_process_id " << (*eventIt)->signal_process_id() ); ATH_MSG_VERBOSE(" signal_process_id " << (*eventIt)->signal_process_id() );
ATH_MSG_VERBOSE(" event number " << (*eventIt)->event_number() ); ATH_MSG_VERBOSE(" event number " << (*eventIt)->event_number() );
ATH_MSG_VERBOSE(" event time " << m_current_event_time << " event index " << m_current_event_index ); ATH_MSG_VERBOSE(" event time " << m_current_event_time << " event index " << m_current_event_index );
...@@ -361,10 +414,10 @@ ISF::GenEventStackFiller::mcEventCollLooper(ISF::ISFParticleContainer& particleC ...@@ -361,10 +414,10 @@ ISF::GenEventStackFiller::mcEventCollLooper(ISF::ISFParticleContainer& particleC
for ( ; passFilter && filterIt!=filterEnd; ++filterIt) { for ( ; passFilter && filterIt!=filterEnd; ++filterIt) {
// determine if the particle passes current filter // determine if the particle passes current filter
passFilter = (*filterIt)->pass(*tParticle); passFilter = (*filterIt)->pass(*tParticle);
ATH_MSG_DEBUG("GenParticleFilter '" << (*filterIt).typeAndName() << "' returned: " ATH_MSG_VERBOSE("GenParticleFilter '" << (*filterIt).typeAndName() << "' returned: "
<< (passFilter ? "true, will keep particle." << (passFilter ? "true, will keep particle."
: "false, will remove particle.")); : "false, will remove particle."));
ATH_MSG_DEBUG("Particle: (" ATH_MSG_VERBOSE("Particle: ("
<<tParticle->momentum().px()<<", " <<tParticle->momentum().px()<<", "
<<tParticle->momentum().py()<<", " <<tParticle->momentum().py()<<", "
<<tParticle->momentum().pz()<<"), pdgCode: " <<tParticle->momentum().pz()<<"), pdgCode: "
...@@ -403,7 +456,7 @@ ISF::GenEventStackFiller::mcEventCollLooper(ISF::ISFParticleContainer& particleC ...@@ -403,7 +456,7 @@ ISF::GenEventStackFiller::mcEventCollLooper(ISF::ISFParticleContainer& particleC
continue; continue;
} }
ATH_MSG_DEBUG( "GenParticle with Barcode '" << pBarcode ATH_MSG_VERBOSE( "GenParticle with Barcode '" << pBarcode
<< "' passed all cuts, adding it to the initial ISF particle list."); << "' passed all cuts, adding it to the initial ISF particle list.");
// -> particle origin (TODO: add proper GeoID, collision/cosmics) // -> particle origin (TODO: add proper GeoID, collision/cosmics)
...@@ -461,10 +514,11 @@ ISF::GenEventStackFiller::mcEventCollLooper(ISF::ISFParticleContainer& particleC ...@@ -461,10 +514,11 @@ ISF::GenEventStackFiller::mcEventCollLooper(ISF::ISFParticleContainer& particleC
sParticle->setUserInformation( new ISF::ParticleUserInformation() ); sParticle->setUserInformation( new ISF::ParticleUserInformation() );
if ( m_barcodeSvc->hasBitCalculator() ) { if ( m_barcodeSvc->hasBitCalculator() ) {
// first gen-event in the list is the hard scatter. // first gen-event in the list is the hard scatter.
if (m_current_event_index==0) { m_barcodeSvc->getBitCalculator()->SetHS(extrabc,true); } // RJH - change thsi to use the pileupCounter parameter
if (nPileupCounter==0) { m_barcodeSvc->getBitCalculator()->SetHS(extrabc,true); }
else { m_barcodeSvc->getBitCalculator()->SetHS(extrabc,false); } else { m_barcodeSvc->getBitCalculator()->SetHS(extrabc,false); }
// bcid // bcid - RJH now set this above from signal-process-id of pileup
int bcid = int( m_current_event_time / m_bunch_spacing ); // int bcid = int( m_current_event_time / m_bunch_spacing );
m_barcodeSvc->getBitCalculator()->SetBCID(extrabc,bcid); m_barcodeSvc->getBitCalculator()->SetBCID(extrabc,bcid);
// parent pid // parent pid
int absPDG = abs( pPdgId ); int absPDG = abs( pPdgId );
...@@ -500,7 +554,7 @@ ISF::GenEventStackFiller::mcEventCollLooper(ISF::ISFParticleContainer& particleC ...@@ -500,7 +554,7 @@ ISF::GenEventStackFiller::mcEventCollLooper(ISF::ISFParticleContainer& particleC
double ISF::GenEventStackFiller::getParticleMass(const HepMC::GenParticle &part) const { double ISF::GenEventStackFiller::getParticleMass(const HepMC::GenParticle &part) const {
// default value: generated particle mass // default value: generated particle mass
double mass = part.generated_mass(); double mass = part.generated_mass();
ATH_MSG_DEBUG("part.generated_mass, mass="<<mass); ATH_MSG_VERBOSE("part.generated_mass, mass="<<mass);
// 1. use PDT mass? // 1. use PDT mass?
if ( !m_useGeneratedParticleMass) { if ( !m_useGeneratedParticleMass) {
...@@ -510,7 +564,7 @@ double ISF::GenEventStackFiller::getParticleMass(const HepMC::GenParticle &part) ...@@ -510,7 +564,7 @@ double ISF::GenEventStackFiller::getParticleMass(const HepMC::GenParticle &part)
: 0 ; : 0 ;
if (pData) { if (pData) {
mass = pData->mass(); mass = pData->mass();
ATH_MSG_DEBUG("using pData mass, mass="<<mass); ATH_MSG_VERBOSE("using pData mass, mass="<<mass);
} }
else else
ATH_MSG_WARNING( "Unable to find mass of particle with PDG ID '" << absPDG << "' in ParticleDataTable. Will set mass to generated_mass: " << mass); ATH_MSG_WARNING( "Unable to find mass of particle with PDG ID '" << absPDG << "' in ParticleDataTable. Will set mass to generated_mass: " << mass);
...@@ -521,6 +575,6 @@ double ISF::GenEventStackFiller::getParticleMass(const HepMC::GenParticle &part) ...@@ -521,6 +575,6 @@ double ISF::GenEventStackFiller::getParticleMass(const HepMC::GenParticle &part)
StatusCode ISF::GenEventStackFiller::finalize() StatusCode ISF::GenEventStackFiller::finalize()
{ {
ATH_MSG_VERBOSE("Finalizing ..."); ATH_MSG_DEBUG("Finalizing ...");
return StatusCode::SUCCESS; return StatusCode::SUCCESS;
} }
...@@ -70,6 +70,7 @@ namespace ISF { ...@@ -70,6 +70,7 @@ namespace ISF {
/** loop over McEventCollection */ /** loop over McEventCollection */
StatusCode processMergedColls(ISF::ISFParticleContainer& particleColl) const; StatusCode processMergedColls(ISF::ISFParticleContainer& particleColl) const;
StatusCode processSingleColl(ISF::ISFParticleContainer& particleColl) const; StatusCode processSingleColl(ISF::ISFParticleContainer& particleColl) const;
StatusCode processPileupColl(ISF::ISFParticleContainer& particleColl) const;
StatusCode mcEventCollLooper(ISF::ISFParticleContainer& particleColl, McEventCollection* mcCollection, int nPileupCounter=0) const; StatusCode mcEventCollLooper(ISF::ISFParticleContainer& particleColl, McEventCollection* mcCollection, int nPileupCounter=0) const;
/** ParticlePropertyService and ParticleDataTable */ /** ParticlePropertyService and ParticleDataTable */
...@@ -78,6 +79,8 @@ namespace ISF { ...@@ -78,6 +79,8 @@ namespace ISF {
std::string m_inputMcEventCollection; //!< name of input McEventCollection std::string m_inputMcEventCollection; //!< name of input McEventCollection
std::string m_outputMcEventCollection; //!< name of output McEventCollection std::string m_outputMcEventCollection; //!< name of output McEventCollection
std::string m_pileupMcEventCollection; //!< name of additional Pileup McEventCollection
std::string m_outputPileupMcEventCollection; //!< name of additional Pileup output McEventCollection
bool m_recordOnlySimulated; //!< record only simulated particles in output truth collection bool m_recordOnlySimulated; //!< record only simulated particles in output truth collection
bool m_useGeneratedParticleMass; //!< use GenParticle::generated_mass() in simulation bool m_useGeneratedParticleMass; //!< use GenParticle::generated_mass() in simulation
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment