9 #include "Pandora/AlgorithmHeaders.h" 27 ThreeDHitCreationAlgorithm::ThreeDHitCreationAlgorithm() :
28 m_iterateTrackHits(true),
29 m_iterateShowerHits(false),
30 m_slidingFitHalfWindow(10),
31 m_nHitRefinementIterations(10),
32 m_sigma3DFitMultiplier(0.2),
33 m_iterationMaxChi2Ratio(1.)
41 for (
const CaloHit *
const pCaloHit : inputCaloHitVector)
43 if (hitType == pCaloHit->GetHitType())
44 outputCaloHitVector.push_back(pCaloHit);
52 const PfoList *pPfoList(
nullptr);
53 PANDORA_RETURN_RESULT_IF_AND_IF(
54 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*
this,
m_inputPfoListName, pPfoList));
56 if (!pPfoList || pPfoList->empty())
58 if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
59 std::cout <<
"ThreeDHitCreationAlgorithm: unable to find pfo list " <<
m_inputPfoListName << std::endl;
61 return STATUS_CODE_SUCCESS;
64 CaloHitList allNewThreeDHits;
66 PfoVector pfoVector(pPfoList->begin(), pPfoList->end());
69 for (
const ParticleFlowObject *
const pPfo : pfoVector)
75 CaloHitVector remainingTwoDHits;
78 if (remainingTwoDHits.empty())
81 pHitCreationTool->Run(
this, pPfo, remainingTwoDHits, protoHitVector);
87 if (protoHitVector.empty())
90 CaloHitList newThreeDHits;
94 allNewThreeDHits.insert(allNewThreeDHits.end(), newThreeDHits.begin(), newThreeDHits.end());
97 if (!allNewThreeDHits.empty())
98 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList(*
this, allNewThreeDHits,
m_outputCaloHitListName));
100 return STATUS_CODE_SUCCESS;
106 const ParticleFlowObject *
const pPfo,
const ProtoHitVector &protoHitVector, CaloHitVector &remainingHitVector)
const 108 ClusterList twoDClusterList;
110 CaloHitList remainingHitList;
112 for (
const Cluster *
const pCluster : twoDClusterList)
115 throw StatusCodeException(STATUS_CODE_FAILURE);
117 pCluster->GetOrderedCaloHitList().FillCaloHitList(remainingHitList);
120 CaloHitSet remainingHitSet(remainingHitList.begin(), remainingHitList.end());
122 for (
const ProtoHit &protoHit : protoHitVector)
126 if (remainingHitSet.end() == eraseIter)
127 throw StatusCodeException(STATUS_CODE_FAILURE);
129 remainingHitSet.erase(eraseIter);
132 remainingHitVector.insert(remainingHitVector.end(), remainingHitSet.begin(), remainingHitSet.end());
143 const float layerPitch(std::max({pitchU, pitchV, pitchW}));
146 double originalChi2(0.);
147 CartesianPointVector currentPoints3D;
148 this->
ExtractResults(protoHitVector, originalChi2, currentPoints3D);
153 const double originalChi2WrtFit(this->
GetChi2WrtFit(originalSlidingFitResult, protoHitVector));
154 double currentChi2(originalChi2 + originalChi2WrtFit);
156 unsigned int nIterations(0);
165 CartesianPointVector newPoints3D;
171 currentChi2 = newChi2;
172 currentPoints3D = newPoints3D;
173 protoHitVector = newProtoHitVector;
176 catch (
const StatusCodeException &)
188 for (
const ProtoHit &protoHit : protoHitVector)
190 chi2 += protoHit.GetChi2();
191 pointVector.push_back(protoHit.GetPosition3D());
202 double chi2WrtFit(0.);
204 for (
const ProtoHit &protoHit : protoHitVector)
206 CartesianVector pointOnFit(0.
f, 0.
f, 0.
f);
212 const double uFit(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoU(pointOnFit.GetY(), pointOnFit.GetZ()));
213 const double vFit(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoV(pointOnFit.GetY(), pointOnFit.GetZ()));
214 const double wFit(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoW(pointOnFit.GetY(), pointOnFit.GetZ()));
216 const double outputU(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoU(
217 protoHit.GetPosition3D().GetY(), protoHit.GetPosition3D().GetZ()));
218 const double outputV(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoV(
219 protoHit.GetPosition3D().GetY(), protoHit.GetPosition3D().GetZ()));
220 const double outputW(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoW(
221 protoHit.GetPosition3D().GetY(), protoHit.GetPosition3D().GetZ()));
223 const double deltaUFit(uFit - outputU), deltaVFit(vFit - outputV), deltaWFit(wFit - outputW);
224 chi2WrtFit += ((deltaUFit * deltaUFit) / (sigma3DFit * sigma3DFit)) + ((deltaVFit * deltaVFit) / (sigma3DFit * sigma3DFit)) +
225 ((deltaWFit * deltaWFit) / (sigma3DFit * sigma3DFit));
236 double hitMovementChi2(0.);
238 for (
const ProtoHit &protoHit : protoHitVector)
240 const CaloHit *
const pCaloHit2D(protoHit.GetParentCaloHit2D());
241 const HitType hitType(pCaloHit2D->GetHitType());
244 const double delta(static_cast<double>(pCaloHit2D->GetPositionVector().GetZ() - projectedPosition.GetZ()));
246 hitMovementChi2 += (delta * delta) / (sigmaUVW * sigmaUVW);
249 return hitMovementChi2;
257 const double sigmaFit(sigmaUVW);
258 const double sigmaHit(sigmaUVW);
261 for (
ProtoHit &protoHit : protoHitVector)
263 CartesianVector pointOnFit(0.
f, 0.
f, 0.
f);
269 const CaloHit *
const pCaloHit2D(protoHit.GetParentCaloHit2D());
270 const HitType hitType(pCaloHit2D->GetHitType());
272 const double uFit(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoU(pointOnFit.GetY(), pointOnFit.GetZ()));
273 const double vFit(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoV(pointOnFit.GetY(), pointOnFit.GetZ()));
274 const double wFit(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoW(pointOnFit.GetY(), pointOnFit.GetZ()));
276 const double sigmaU((TPC_VIEW_U == hitType) ? sigmaHit : sigmaFit);
277 const double sigmaV((TPC_VIEW_V == hitType) ? sigmaHit : sigmaFit);
278 const double sigmaW((TPC_VIEW_W == hitType) ? sigmaHit : sigmaFit);
280 CartesianVector position3D(0.
f, 0.
f, 0.
f);
281 double chi2(std::numeric_limits<double>::max());
282 double u(std::numeric_limits<double>::max()), v(std::numeric_limits<double>::max()),
w(std::numeric_limits<double>::max());
284 if (protoHit.GetNTrajectorySamples() == 2)
286 u = (TPC_VIEW_U == hitType) ? pCaloHit2D->GetPositionVector().GetZ()
287 : (TPC_VIEW_U == protoHit.GetFirstTrajectorySample().GetHitType()) ? protoHit.GetFirstTrajectorySample().GetPosition().GetZ()
288 : protoHit.GetLastTrajectorySample().GetPosition().GetZ();
289 v = (TPC_VIEW_V == hitType) ? pCaloHit2D->GetPositionVector().GetZ()
290 : (TPC_VIEW_V == protoHit.GetFirstTrajectorySample().GetHitType()) ? protoHit.GetFirstTrajectorySample().GetPosition().GetZ()
291 : protoHit.GetLastTrajectorySample().GetPosition().GetZ();
292 w = (TPC_VIEW_W == hitType) ? pCaloHit2D->GetPositionVector().GetZ()
293 : (TPC_VIEW_W == protoHit.GetFirstTrajectorySample().GetHitType()) ? protoHit.GetFirstTrajectorySample().GetPosition().GetZ()
294 : protoHit.GetLastTrajectorySample().GetPosition().GetZ();
296 else if (protoHit.GetNTrajectorySamples() == 1)
298 u = PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoU(
299 protoHit.GetPosition3D().GetY(), protoHit.GetPosition3D().GetZ());
300 v = PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoV(
301 protoHit.GetPosition3D().GetY(), protoHit.GetPosition3D().GetZ());
302 w = PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoW(
303 protoHit.GetPosition3D().GetY(), protoHit.GetPosition3D().GetZ());
307 std::cout <<
"ThreeDHitCreationAlgorithm::IterativeTreatment - Unexpected number of trajectory samples" << std::endl;
308 throw StatusCodeException(STATUS_CODE_FAILURE);
311 double bestY(std::numeric_limits<double>::max()), bestZ(std::numeric_limits<double>::max());
312 PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->GetMinChiSquaredYZ(
313 u, v, w, sigmaU, sigmaV, sigmaW, uFit, vFit, wFit, sigma3DFit, bestY, bestZ, chi2);
314 position3D.SetValues(pCaloHit2D->GetPositionVector().GetX(),
static_cast<float>(bestY), static_cast<float>(bestZ));
316 protoHit.SetPosition3D(position3D, chi2);
324 for (
const ProtoHit &protoHit : protoHitVector)
326 const CaloHit *pCaloHit3D(
nullptr);
330 throw StatusCodeException(STATUS_CODE_FAILURE);
332 newThreeDHits.push_back(pCaloHit3D);
341 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
343 PandoraContentApi::CaloHit::Parameters parameters;
345 parameters.m_hitType = TPC_3D;
348 parameters.m_pParentAddress =
static_cast<const void *
>(pCaloHit2D);
351 parameters.m_cellThickness = pCaloHit2D->GetCellThickness();
352 parameters.m_cellGeometry = RECTANGULAR;
353 parameters.m_cellSize0 = pCaloHit2D->GetCellLengthScale();
354 parameters.m_cellSize1 = pCaloHit2D->GetCellLengthScale();
355 parameters.m_cellNormalVector = pCaloHit2D->GetCellNormalVector();
356 parameters.m_expectedDirection = pCaloHit2D->GetExpectedDirection();
357 parameters.m_nCellRadiationLengths = pCaloHit2D->GetNCellRadiationLengths();
358 parameters.m_nCellInteractionLengths = pCaloHit2D->GetNCellInteractionLengths();
359 parameters.m_time = pCaloHit2D->GetTime();
360 parameters.m_inputEnergy = pCaloHit2D->GetInputEnergy();
361 parameters.m_mipEquivalentEnergy = pCaloHit2D->GetMipEquivalentEnergy();
362 parameters.m_electromagneticEnergy = pCaloHit2D->GetElectromagneticEnergy();
363 parameters.m_hadronicEnergy = pCaloHit2D->GetHadronicEnergy();
364 parameters.m_isDigital = pCaloHit2D->IsDigital();
365 parameters.m_hitRegion = pCaloHit2D->GetHitRegion();
366 parameters.m_layer = pCaloHit2D->GetLayer();
367 parameters.m_isInOuterSamplingLayer = pCaloHit2D->IsInOuterSamplingLayer();
368 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CaloHit::Create(*
this, parameters, pCaloHit3D));
378 (void)PandoraContentApi::GetPlugins(*this)->GetPseudoLayerPlugin()->GetPseudoLayer(protoHit.
GetPosition3D());
380 catch (StatusCodeException &)
393 if (caloHitList.empty())
394 throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
396 ClusterList threeDClusterList;
399 if (!threeDClusterList.empty())
400 throw StatusCodeException(STATUS_CODE_FAILURE);
402 const ClusterList *pClusterList(
nullptr);
403 std::string clusterListName;
404 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*
this, pClusterList, clusterListName));
406 PandoraContentApi::Cluster::Parameters parameters;
407 parameters.m_caloHitList.insert(parameters.m_caloHitList.end(), caloHitList.begin(), caloHitList.end());
409 const Cluster *pCluster3D(
nullptr);
410 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Cluster::Create(*
this, parameters, pCluster3D));
412 if (!pCluster3D || !pClusterList || pClusterList->empty())
413 throw StatusCodeException(STATUS_CODE_FAILURE);
415 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Cluster>(*
this,
m_outputClusterListName));
416 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToPfo(*
this, pPfo, pCluster3D));
424 if (!m_isPositionSet)
425 throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
434 if (!m_isPositionSet)
435 throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
444 if (m_trajectorySampleVector.empty())
445 throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
447 return m_trajectorySampleVector.front();
454 if (m_trajectorySampleVector.size() < 2)
455 throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
457 return m_trajectorySampleVector.back();
465 AlgorithmToolVector algorithmToolVector;
466 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ProcessAlgorithmToolList(*
this, xmlHandle,
"HitCreationTools", algorithmToolVector));
472 if (!pHitCreationTool)
473 return STATUS_CODE_INVALID_PARAMETER;
478 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"InputPfoListName",
m_inputPfoListName));
479 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"OutputCaloHitListName",
m_outputCaloHitListName));
480 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"OutputClusterListName",
m_outputClusterListName));
482 PANDORA_RETURN_RESULT_IF_AND_IF(
483 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"IterateTrackHits",
m_iterateTrackHits));
485 PANDORA_RETURN_RESULT_IF_AND_IF(
486 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"IterateShowerHits",
m_iterateShowerHits));
488 PANDORA_RETURN_RESULT_IF_AND_IF(
489 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"SlidingFitHalfWindow",
m_slidingFitHalfWindow));
491 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
494 PANDORA_RETURN_RESULT_IF_AND_IF(
495 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"Sigma3DFitMultiplier",
m_sigma3DFitMultiplier));
497 PANDORA_RETURN_RESULT_IF_AND_IF(
498 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"IterationMaxChi2Ratio",
m_iterationMaxChi2Ratio));
500 return STATUS_CODE_SUCCESS;
void IterativeTreatment(ProtoHitVector &protoHitVector) const
Improve initial 3D hits by fitting proto hits and iteratively creating consisted 3D hit trajectory...
static bool SortByNHits(const pandora::ParticleFlowObject *const pLhs, const pandora::ParticleFlowObject *const pRhs)
Sort pfos by number of constituent hits.
Header file for the pfo helper class.
Proto hits are temporary constructs to be used during iterative 3D hit procedure. ...
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
static float GetSigmaUVW(const pandora::Pandora &pandora, const float maxSigmaDiscrepancy=0.01)
Find the sigmaUVW value for the detector geometry.
std::string m_inputPfoListName
The name of the input pfo list.
unsigned int m_slidingFitHalfWindow
The sliding linear fit half window.
bool CheckThreeDHit(const ProtoHit &protoHit) const
Check that a new three dimensional position is not unphysical.
static void GetTwoDClusterList(const pandora::ParticleFlowObject *const pPfo, pandora::ClusterList &clusterList)
Get the list of 2D clusters from an input pfo.
const pandora::CaloHit * GetParentCaloHit2D() const
Get the address of the parent 2D calo hit.
std::string m_outputCaloHitListName
The name of the output calo hit list.
void SeparateTwoDHits(const pandora::ParticleFlowObject *const pPfo, const ProtoHitVector &protoHitVector, pandora::CaloHitVector &remainingHitVector) const
Get the list of 2D calo hits in a pfo for which 3D hits have and have not been created.
static pandora::CartesianVector ProjectPosition(const pandora::Pandora &pandora, const pandora::CartesianVector &position3D, const pandora::HitType view)
Project 3D position into a given 2D view.
void RefineHitPositions(const ThreeDSlidingFitResult &slidingFitResult, ProtoHitVector &protoHitVector) const
Refine the 3D hit positions (and chi2) for a list of proto hits, in accordance with a provided 3D sli...
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
void CreateThreeDHit(const ProtoHit &protoHit, const pandora::CaloHit *&pCaloHit3D) const
Create a new three dimensional hit from a two dimensional hit.
static bool IsTrack(const pandora::ParticleFlowObject *const pPfo)
Return track flag based on Pfo Particle ID.
double GetHitMovementChi2(const ProtoHitVector &protoHitVector) const
Receive a chi2 value indicating consistency of a list of proto hits with the original, input hit positions.
float GetLongitudinalDisplacement(const pandora::CartesianVector &position) const
Get longitudinal projection onto primary axis.
unsigned int m_nHitRefinementIterations
The maximum number of hit refinement iterations.
const pandora::CartesianVector & GetPosition3D() const
Get the output 3D position.
Header file for the three dimensional hit creation algorithm class.
Header file for the geometry helper class.
std::string m_outputClusterListName
The name of the output cluster list.
pandora::StatusCode Run()
static bool IsShower(const pandora::ParticleFlowObject *const pPfo)
Return shower flag based on Pfo Particle ID.
double GetChi2WrtFit(const ThreeDSlidingFitResult &slidingFitResult, const ProtoHitVector &protoHitVector) const
Receive a chi2 value indicating consistency of a list of proto hits with a provided 3D sliding fit tr...
double GetChi2() const
Get the chi squared value.
static bool SortHitsByPosition(const pandora::CaloHit *const pLhs, const pandora::CaloHit *const pRhs)
Sort calo hits by their position (use Z, followed by X, followed by Y)
Header file for the cluster helper class.
bool m_iterateTrackHits
Whether to enable iterative improvement of 3D hits for track trajectories.
std::vector< ProtoHit > ProtoHitVector
bool m_iterateShowerHits
Whether to enable iterative improvement of 3D hits for showers.
static float GetWirePitch(const pandora::Pandora &pandora, const pandora::HitType view, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
const TrajectorySample & GetFirstTrajectorySample() const
Get the first trajectory sample.
static void GetThreeDClusterList(const pandora::ParticleFlowObject *const pPfo, pandora::ClusterList &clusterList)
Get the list of 3D clusters from an input pfo.
void CreateThreeDHits(const ProtoHitVector &protoHitVector, pandora::CaloHitList &newThreeDHits) const
Create new three dimensional hits from two dimensional hits.
Header file for the lar three dimensional sliding fit result class.
ThreeDSlidingFitResult class.
HitCreationBaseTool class.
void ExtractResults(const ProtoHitVector &protoHitVector, double &chi2, pandora::CartesianPointVector &pointVector) const
Extract key results from a provided proto hit vector.
Trajectory samples record the results of sampling a particles in a particular view.
pandora::StatusCode GetGlobalFitPosition(const float rL, pandora::CartesianVector &position) const
Get global fit position for a given longitudinal coordinate.
HitCreationToolVector m_algorithmToolVector
The algorithm tool vector.
double m_sigma3DFitMultiplier
Multiplicative factor: sigmaUVW (same as sigmaHit and sigma2DFit) to sigma3DFit.
double m_iterationMaxChi2Ratio
Max ratio between current and previous chi2 values to cease iterations.
void AddThreeDHitsToPfo(const pandora::ParticleFlowObject *const pPfo, const pandora::CaloHitList &caloHitList) const
Add a specified list of three dimensional hits to a cluster in a pfo, creating the new cluster if req...
const TrajectorySample & GetLastTrajectorySample() const
Get the last trajectory sample.
void FilterCaloHitsByType(const pandora::CaloHitVector &inputCaloHitVector, const pandora::HitType hitType, pandora::CaloHitVector &outputCaloHitVector) const
Get the subset of a provided calo hit vector corresponding to a specified hit type.