diff --git a/source/Utils/include/FilterConeHits.h b/source/Utils/include/FilterConeHits.h index c1b31a6..9c968e1 100644 --- a/source/Utils/include/FilterConeHits.h +++ b/source/Utils/include/FilterConeHits.h @@ -24,6 +24,7 @@ using namespace marlin; * @parameter TrackerHitOutputRelations name of the tracker hit relation output collections * @parameter DeltaRCut maximum angular distance between the hits and the particle direction * @parameter FillHistograms flag to fill the diagnostic histograms + * @parameter ConeAroundStatus list of generator statuses that we can choose to cone around * * @author M. Casarsa, INFN Trieste * @date 22 January 2021 @@ -63,6 +64,9 @@ class FilterConeHits : public Processor { double m_deltaRCut{}; double m_dist3DCut{}; + // Coning around MCParticles whose generator status is in this list, given in steering file + std::vector m_coneAroundStatus{}; + // --- Diagnostic histograms: TH1F* m_distXY = nullptr; TH1F* m_distZ = nullptr; diff --git a/source/Utils/src/FilterConeHits.cc b/source/Utils/src/FilterConeHits.cc index 6c79c55..cf6339a 100644 --- a/source/Utils/src/FilterConeHits.cc +++ b/source/Utils/src/FilterConeHits.cc @@ -2,6 +2,7 @@ #include #include #include +#include #include @@ -29,7 +30,6 @@ FilterConeHits::FilterConeHits() : Processor("FilterConeHits") { _description = "FilterConeHits selects tracker hits in a cone opened around a MC particle direction"; // --- Processor parameters: - registerProcessorParameter("MCParticleCollection", "Name of the MCParticle collection", m_inputMCParticlesCollName, std::string("MCParticle")); @@ -58,6 +58,12 @@ FilterConeHits::FilterConeHits() : Processor("FilterConeHits") { double(-1.)); registerProcessorParameter("FillHistograms", "Flag to fill the diagnostic histograms", m_fillHistos, false); + + registerProcessorParameter( "ConeAroundStatus", + "List of MC Particle statuses to build cones around. Default is 1", + m_coneAroundStatus, + std::vector{1} ); + } void FilterConeHits::init() { @@ -91,6 +97,8 @@ void FilterConeHits::init() { void FilterConeHits::processRunHeader(LCRunHeader*) { _nRun++; } void FilterConeHits::processEvent(LCEvent* evt) { + + streamlog_out(DEBUG) << "[FilterConeHits] processEvent run=" << evt->getRunNumber() << " evt=" << evt->getEventNumber() << std::endl; // --- Check whether the number of input and output collections match if (m_inputTrackerHitsCollNames.size() != m_inputTrackerSimHitsCollNames.size() || @@ -146,6 +154,7 @@ void FilterConeHits::processEvent(LCEvent* evt) { continue; } + // reco hit output collections std::string encoderString = inputHitColls[icol]->getParameters().getStringVal("CellIDEncoding"); outputTrackerHitColls[icol] = new LCCollectionVec(inputHitColls[icol]->getTypeName()); @@ -181,13 +190,14 @@ void FilterConeHits::processEvent(LCEvent* evt) { for (int ipart = 0; ipart < m_inputMCParticles->getNumberOfElements(); ++ipart) { MCParticle* part = dynamic_cast(m_inputMCParticles->getElementAt(ipart)); - // --- Keep only the generator-level particles: - if (part->getGeneratorStatus() != 1) + const int genStat = part->getGeneratorStatus(); + if ( std::find(m_coneAroundStatus.begin(), m_coneAroundStatus.end(), genStat) == m_coneAroundStatus.end() ) { continue; - - double part_p = - sqrt(part->getMomentum()[0] * part->getMomentum()[0] + part->getMomentum()[1] * part->getMomentum()[1] + - part->getMomentum()[2] * part->getMomentum()[2]); + } + + double part_p = sqrt( part->getMomentum()[0]*part->getMomentum()[0] + + part->getMomentum()[1]*part->getMomentum()[1] + + part->getMomentum()[2]*part->getMomentum()[2] ); HelixClass_double helix; helix.Initialize_VP((double*)part->getVertex(), (double*)part->getMomentum(), (double)part->getCharge(), @@ -280,7 +290,6 @@ void FilterConeHits::processEvent(LCEvent* evt) { << outputTrackerHitColls[icol]->getTypeName() << " added to the event \n" << " output collection " << m_outputTrackerSimHitsCollNames[icol] << " of type " << outputTrackerSimHitColls[icol]->getTypeName() << " added to the event " << std::endl; - } // icol loop streamlog_out(DEBUG) << " processing event: " << evt->getEventNumber() << " in run: " << evt->getRunNumber()