HitCreationBaseTool.cc

Go to the documentation of this file.

 
#include "Pandora/AlgorithmHeaders.h"
 
#include "larpandoracontent/LArHelpers/LArGeometryHelper.h"
 
#include "larpandoracontent/LArThreeDReco/LArHitCreation/HitCreationBaseTool.h"
 
using namespace pandora;
 
namespace lar_content
{
 
HitCreationBaseTool::HitCreationBaseTool() : m_sigmaX2(1.), m_chiSquaredCut(1.)
{
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
HitCreationBaseTool::~HitCreationBaseTool()
{
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void HitCreationBaseTool::GetBestPosition3D(const HitType hitType1, const HitType hitType2, const CartesianPointVector &fitPositionList1,
    const CartesianPointVector &fitPositionList2, ProtoHit &protoHit) const
{
    if (fitPositionList1.empty() && fitPositionList2.empty())
    {
        throw StatusCodeException(STATUS_CODE_NOT_FOUND);
    }
    else if (fitPositionList1.empty())
    {
        if (fitPositionList2.size() != 1)
            throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
 
        this->GetBestPosition3D(hitType2, fitPositionList2.front(), protoHit);
    }
    else if (fitPositionList2.empty())
    {
        if (fitPositionList1.size() != 1)
            throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
 
        this->GetBestPosition3D(hitType1, fitPositionList1.front(), protoHit);
    }
    else
    {
        for (const CartesianVector &fitPosition1 : fitPositionList1)
        {
            for (const CartesianVector &fitPosition2 : fitPositionList2)
            {
                ProtoHit thisProtoHit(protoHit.GetParentCaloHit2D());
                this->GetBestPosition3D(hitType1, hitType2, fitPosition1, fitPosition2, thisProtoHit);
 
                if (!protoHit.IsPositionSet() || (thisProtoHit.GetChi2() < protoHit.GetChi2()))
                    protoHit = thisProtoHit;
            }
        }
    }
 
    // TODO Add additional chi-squared at edge of detector
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void HitCreationBaseTool::GetBestPosition3D(const HitType hitType1, const HitType hitType2, const CartesianVector &fitPosition1,
    const CartesianVector &fitPosition2, ProtoHit &protoHit) const
{
    // TODO Input better uncertainties into this method (sigmaHit, sigmaFit, sigmaX)
    const CaloHit *const pCaloHit2D(protoHit.GetParentCaloHit2D());
    const HitType hitType(pCaloHit2D->GetHitType());
 
    const double sigmaFit(LArGeometryHelper::GetSigmaUVW(this->GetPandora()));
    const double sigmaHit(sigmaFit);
 
    CartesianVector position3D(0.f, 0.f, 0.f);
    double chi2(std::numeric_limits<double>::max());
 
    const double u((TPC_VIEW_U == hitType) ? pCaloHit2D->GetPositionVector().GetZ()
                                           : (TPC_VIEW_U == hitType1) ? fitPosition1.GetZ() : fitPosition2.GetZ());
    const double v((TPC_VIEW_V == hitType) ? pCaloHit2D->GetPositionVector().GetZ()
                                           : (TPC_VIEW_V == hitType1) ? fitPosition1.GetZ() : fitPosition2.GetZ());
    const double w((TPC_VIEW_W == hitType) ? pCaloHit2D->GetPositionVector().GetZ()
                                           : (TPC_VIEW_W == hitType1) ? fitPosition1.GetZ() : fitPosition2.GetZ());
 
    const double sigmaU((TPC_VIEW_U == hitType) ? sigmaHit : sigmaFit);
    const double sigmaV((TPC_VIEW_V == hitType) ? sigmaHit : sigmaFit);
    const double sigmaW((TPC_VIEW_W == hitType) ? sigmaHit : sigmaFit);
 
    double bestY(std::numeric_limits<double>::max()), bestZ(std::numeric_limits<double>::max());
    this->GetPandora().GetPlugins()->GetLArTransformationPlugin()->GetMinChiSquaredYZ(u, v, w, sigmaU, sigmaV, sigmaW, bestY, bestZ, chi2);
    position3D.SetValues(pCaloHit2D->GetPositionVector().GetX(), static_cast<float>(bestY), static_cast<float>(bestZ));
 
    const double deltaX1(pCaloHit2D->GetPositionVector().GetX() - fitPosition1.GetX());
    const double deltaX2(pCaloHit2D->GetPositionVector().GetX() - fitPosition2.GetX());
    const double chi2X(((deltaX1 * deltaX1) / m_sigmaX2) + ((deltaX2 * deltaX2) / m_sigmaX2));
 
    protoHit.SetPosition3D(position3D, chi2 + chi2X);
    protoHit.AddTrajectorySample(TrajectorySample(fitPosition1, hitType1, sigmaFit));
    protoHit.AddTrajectorySample(TrajectorySample(fitPosition2, hitType2, sigmaFit));
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void HitCreationBaseTool::GetBestPosition3D(const HitType hitType, const CartesianVector &fitPosition, ProtoHit &protoHit) const
{
    // TODO Input better uncertainties into this method (sigmaHit, sigmaFit, sigmaX)
    const CaloHit *const pCaloHit2D(protoHit.GetParentCaloHit2D());
    const double sigmaFit(LArGeometryHelper::GetSigmaUVW(this->GetPandora()));
 
    if (pCaloHit2D->GetHitType() == hitType)
        throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
 
    CartesianVector position3D(0.f, 0.f, 0.f);
    float chi2(std::numeric_limits<float>::max());
    LArGeometryHelper::MergeTwoPositions3D(
        this->GetPandora(), pCaloHit2D->GetHitType(), hitType, pCaloHit2D->GetPositionVector(), fitPosition, position3D, chi2);
 
    // ATTN Replace chi2 from LArGeometryHelper for consistency with three-view treatment (purely a measure of delta x)
    const double deltaX(pCaloHit2D->GetPositionVector().GetX() - fitPosition.GetX());
    const double chi2X((deltaX * deltaX) / m_sigmaX2);
 
    protoHit.SetPosition3D(position3D, chi2X);
    protoHit.AddTrajectorySample(TrajectorySample(fitPosition, hitType, sigmaFit));
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
StatusCode HitCreationBaseTool::ReadSettings(const pandora::TiXmlHandle xmlHandle)
{
    double sigmaX(std::sqrt(m_sigmaX2));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SigmaX", sigmaX));
 
    m_sigmaX2 = sigmaX * sigmaX;
 
    if (m_sigmaX2 < std::numeric_limits<double>::epsilon())
    {
        std::cout << "HitCreationBaseTool - Invalid parameter, SigmaX: " << sigmaX << std::endl;
        return STATUS_CODE_INVALID_PARAMETER;
    }
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ChiSquaredCut", m_chiSquaredCut));
 
    return STATUS_CODE_SUCCESS;
}
 
} // namespace lar_content