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(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*
this,
m_inputPfoListName, pPfoList));
55 if (!pPfoList || pPfoList->empty())
57 if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
58 std::cout <<
"ThreeDHitCreationAlgorithm: unable to find pfo list " <<
m_inputPfoListName << std::endl;
60 return STATUS_CODE_SUCCESS;
63 CaloHitList allNewThreeDHits;
65 PfoVector pfoVector(pPfoList->begin(), pPfoList->end());
68 for (
const ParticleFlowObject *
const pPfo : pfoVector)
74 CaloHitVector remainingTwoDHits;
77 if (remainingTwoDHits.empty())
80 pHitCreationTool->Run(
this, pPfo, remainingTwoDHits, protoHitVector);
86 if (protoHitVector.empty())
89 CaloHitList newThreeDHits;
93 allNewThreeDHits.insert(allNewThreeDHits.end(), newThreeDHits.begin(), newThreeDHits.end());
96 if (!allNewThreeDHits.empty())
97 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList(*
this, allNewThreeDHits,
m_outputCaloHitListName));
99 return STATUS_CODE_SUCCESS;
106 ClusterList twoDClusterList;
108 CaloHitList remainingHitList;
110 for (
const Cluster *
const pCluster : twoDClusterList)
113 throw StatusCodeException(STATUS_CODE_FAILURE);
115 pCluster->GetOrderedCaloHitList().FillCaloHitList(remainingHitList);
118 CaloHitSet remainingHitSet(remainingHitList.begin(), remainingHitList.end());
120 for (
const ProtoHit &protoHit : protoHitVector)
124 if (remainingHitSet.end() == eraseIter)
125 throw StatusCodeException(STATUS_CODE_FAILURE);
127 remainingHitSet.erase(eraseIter);
130 remainingHitVector.insert(remainingHitVector.end(), remainingHitSet.begin(), remainingHitSet.end());
141 double originalChi2(0.);
142 CartesianPointVector currentPoints3D;
143 this->
ExtractResults(protoHitVector, originalChi2, currentPoints3D);
148 const double originalChi2WrtFit(this->
GetChi2WrtFit(originalSlidingFitResult, protoHitVector));
149 double currentChi2(originalChi2 + originalChi2WrtFit);
151 unsigned int nIterations(0);
160 CartesianPointVector newPoints3D;
166 currentChi2 = newChi2;
167 currentPoints3D = newPoints3D;
168 protoHitVector = newProtoHitVector;
171 catch (
const StatusCodeException &)
183 for (
const ProtoHit &protoHit : protoHitVector)
185 chi2 += protoHit.GetChi2();
186 pointVector.push_back(protoHit.GetPosition3D());
197 double chi2WrtFit(0.);
199 for (
const ProtoHit &protoHit : protoHitVector)
201 CartesianVector pointOnFit(0.
f, 0.
f, 0.
f);
207 const double uFit(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoU(pointOnFit.GetY(), pointOnFit.GetZ()));
208 const double vFit(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoV(pointOnFit.GetY(), pointOnFit.GetZ()));
209 const double wFit(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoW(pointOnFit.GetY(), pointOnFit.GetZ()));
211 const double outputU(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoU(protoHit.GetPosition3D().GetY(), protoHit.GetPosition3D().GetZ()));
212 const double outputV(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoV(protoHit.GetPosition3D().GetY(), protoHit.GetPosition3D().GetZ()));
213 const double outputW(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoW(protoHit.GetPosition3D().GetY(), protoHit.GetPosition3D().GetZ()));
215 const double deltaUFit(uFit - outputU), deltaVFit(vFit - outputV), deltaWFit(wFit - outputW);
216 chi2WrtFit += ((deltaUFit * deltaUFit) / (sigma3DFit * sigma3DFit)) + ((deltaVFit * deltaVFit) / (sigma3DFit * sigma3DFit)) + ((deltaWFit * deltaWFit) / (sigma3DFit * sigma3DFit));
227 double hitMovementChi2(0.);
229 for (
const ProtoHit &protoHit : protoHitVector)
231 const CaloHit *
const pCaloHit2D(protoHit.GetParentCaloHit2D());
232 const HitType hitType(pCaloHit2D->GetHitType());
235 const double delta(static_cast<double>(pCaloHit2D->GetPositionVector().GetZ() - projectedPosition.GetZ()));
237 hitMovementChi2 += (delta * delta) / (sigmaUVW * sigmaUVW);
240 return hitMovementChi2;
248 const double sigmaFit(sigmaUVW);
249 const double sigmaHit(sigmaUVW);
252 for (
ProtoHit &protoHit : protoHitVector)
254 CartesianVector pointOnFit(0.
f, 0.
f, 0.
f);
260 const CaloHit *
const pCaloHit2D(protoHit.GetParentCaloHit2D());
261 const HitType hitType(pCaloHit2D->GetHitType());
263 const double uFit(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoU(pointOnFit.GetY(), pointOnFit.GetZ()));
264 const double vFit(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoV(pointOnFit.GetY(), pointOnFit.GetZ()));
265 const double wFit(PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoW(pointOnFit.GetY(), pointOnFit.GetZ()));
267 const double sigmaU((TPC_VIEW_U == hitType) ? sigmaHit : sigmaFit);
268 const double sigmaV((TPC_VIEW_V == hitType) ? sigmaHit : sigmaFit);
269 const double sigmaW((TPC_VIEW_W == hitType) ? sigmaHit : sigmaFit);
271 CartesianVector position3D(0.
f, 0.
f, 0.
f);
275 if (protoHit.GetNTrajectorySamples() == 2)
277 u = (TPC_VIEW_U == hitType) ? pCaloHit2D->GetPositionVector().GetZ() : (TPC_VIEW_U == protoHit.GetFirstTrajectorySample().GetHitType()) ? protoHit.GetFirstTrajectorySample().GetPosition().GetZ() : protoHit.GetLastTrajectorySample().GetPosition().GetZ();
278 v = (TPC_VIEW_V == hitType) ? pCaloHit2D->GetPositionVector().GetZ() : (TPC_VIEW_V == protoHit.GetFirstTrajectorySample().GetHitType()) ? protoHit.GetFirstTrajectorySample().GetPosition().GetZ() : protoHit.GetLastTrajectorySample().GetPosition().GetZ();
279 w = (TPC_VIEW_W == hitType) ? pCaloHit2D->GetPositionVector().GetZ() : (TPC_VIEW_W == protoHit.GetFirstTrajectorySample().GetHitType()) ? protoHit.GetFirstTrajectorySample().GetPosition().GetZ() : protoHit.GetLastTrajectorySample().GetPosition().GetZ();
281 else if (protoHit.GetNTrajectorySamples() == 1)
283 u = PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoU(protoHit.GetPosition3D().GetY(), protoHit.GetPosition3D().GetZ());
284 v = PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoV(protoHit.GetPosition3D().GetY(), protoHit.GetPosition3D().GetZ());
285 w = PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->YZtoW(protoHit.GetPosition3D().GetY(), protoHit.GetPosition3D().GetZ());
289 std::cout <<
"ThreeDHitCreationAlgorithm::IterativeTreatment - Unexpected number of trajectory samples" << std::endl;
290 throw StatusCodeException(STATUS_CODE_FAILURE);
294 PandoraContentApi::GetPlugins(*this)->GetLArTransformationPlugin()->GetMinChiSquaredYZ(u, v, w, sigmaU, sigmaV, sigmaW, uFit, vFit, wFit, sigma3DFit, bestY, bestZ, chi2);
295 position3D.SetValues(pCaloHit2D->GetPositionVector().GetX(),
static_cast<float>(bestY), static_cast<float>(bestZ));
297 protoHit.SetPosition3D(position3D, chi2);
305 for (
const ProtoHit &protoHit : protoHitVector)
307 const CaloHit *pCaloHit3D(
nullptr);
311 throw StatusCodeException(STATUS_CODE_FAILURE);
313 newThreeDHits.push_back(pCaloHit3D);
322 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
324 PandoraContentApi::CaloHit::Parameters parameters;
326 parameters.m_hitType = TPC_3D;
329 parameters.m_pParentAddress =
static_cast<const void*
>(pCaloHit2D);
332 parameters.m_cellThickness = pCaloHit2D->GetCellThickness();
333 parameters.m_cellGeometry = RECTANGULAR;
334 parameters.m_cellSize0 = pCaloHit2D->GetCellLengthScale();
335 parameters.m_cellSize1 = pCaloHit2D->GetCellLengthScale();
336 parameters.m_cellNormalVector = pCaloHit2D->GetCellNormalVector();
337 parameters.m_expectedDirection = pCaloHit2D->GetExpectedDirection();
338 parameters.m_nCellRadiationLengths = pCaloHit2D->GetNCellRadiationLengths();
339 parameters.m_nCellInteractionLengths = pCaloHit2D->GetNCellInteractionLengths();
340 parameters.m_time = pCaloHit2D->GetTime();
341 parameters.m_inputEnergy = pCaloHit2D->GetInputEnergy();
342 parameters.m_mipEquivalentEnergy = pCaloHit2D->GetMipEquivalentEnergy();
343 parameters.m_electromagneticEnergy = pCaloHit2D->GetElectromagneticEnergy();
344 parameters.m_hadronicEnergy = pCaloHit2D->GetHadronicEnergy();
345 parameters.m_isDigital = pCaloHit2D->IsDigital();
346 parameters.m_hitRegion = pCaloHit2D->GetHitRegion();
347 parameters.m_layer = pCaloHit2D->GetLayer();
348 parameters.m_isInOuterSamplingLayer = pCaloHit2D->IsInOuterSamplingLayer();
349 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CaloHit::Create(*
this, parameters, pCaloHit3D));
359 (void) PandoraContentApi::GetPlugins(*this)->GetPseudoLayerPlugin()->GetPseudoLayer(protoHit.
GetPosition3D());
361 catch (StatusCodeException &)
374 if (caloHitList.empty())
375 throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
377 ClusterList threeDClusterList;
380 if (!threeDClusterList.empty())
381 throw StatusCodeException(STATUS_CODE_FAILURE);
383 const ClusterList *pClusterList(
nullptr); std::string clusterListName;
384 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*
this, pClusterList, clusterListName));
386 PandoraContentApi::Cluster::Parameters parameters;
387 parameters.m_caloHitList.insert(parameters.m_caloHitList.end(), caloHitList.begin(), caloHitList.end());
389 const Cluster *pCluster3D(
nullptr);
390 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Cluster::Create(*
this, parameters, pCluster3D));
392 if (!pCluster3D || !pClusterList || pClusterList->empty())
393 throw StatusCodeException(STATUS_CODE_FAILURE);
395 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Cluster>(*
this,
m_outputClusterListName));
396 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToPfo(*
this, pPfo, pCluster3D));
404 if (!m_isPositionSet)
405 throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
414 if (!m_isPositionSet)
415 throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
424 if (m_trajectorySampleVector.empty())
425 throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
427 return m_trajectorySampleVector.front();
434 if (m_trajectorySampleVector.size() < 2)
435 throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
437 return m_trajectorySampleVector.back();
445 AlgorithmToolVector algorithmToolVector;
446 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ProcessAlgorithmToolList(*
this, xmlHandle,
447 "HitCreationTools", algorithmToolVector));
453 if (!pHitCreationTool)
454 return STATUS_CODE_INVALID_PARAMETER;
459 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"InputPfoListName",
m_inputPfoListName));
460 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"OutputCaloHitListName",
m_outputCaloHitListName));
461 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"OutputClusterListName",
m_outputClusterListName));
463 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
466 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
469 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
472 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
475 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
478 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
481 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.
static float GetWireZPitch(const pandora::Pandora &pandora, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
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.
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.