EnergyDepositionAsymmetryFeatureTool.cc

Go to the documentation of this file.

#include "larpandoracontent/LArVertex/EnergyDepositionAsymmetryFeatureTool.h"
#include "Pandora/AlgorithmHeaders.h"
#include "larpandoracontent/LArHelpers/LArClusterHelper.h"
#include "larpandoracontent/LArHelpers/LArGeometryHelper.h"

using namespace pandora;

namespace lar_content
{

EnergyDepositionAsymmetryFeatureTool::EnergyDepositionAsymmetryFeatureTool() :
    GlobalAsymmetryFeatureTool()
{
}

//------------------------------------------------------------------------------------------------------------------------------------------

float EnergyDepositionAsymmetryFeatureTool::CalculateAsymmetry(const bool useEnergyMetrics, const CartesianVector &vertexPosition2D,
    const ClusterVector &clusterVector, const CartesianVector &localWeightedDirectionSum) const
{
    // Project every hit onto local event axis direction and record side of the projected vtx position on which it falls
    float beforeVtxEnergy(0.f), afterVtxEnergy(0.f);
    unsigned int beforeVtxHits(0), afterVtxHits(0);

    const CartesianVector localWeightedDirection(localWeightedDirectionSum.GetUnitVector());
    const float evtProjectedVtxPos(vertexPosition2D.GetDotProduct(localWeightedDirection));

    float minBeforeProjectedPos(std::numeric_limits<float>::max());
    float maxBeforeProjectedPos(-std::numeric_limits<float>::max());

    float minAfterProjectedPos(std::numeric_limits<float>::max());
    float maxAfterProjectedPos(-std::numeric_limits<float>::max());

    for (const Cluster *const pCluster : clusterVector)
    {
        CaloHitList caloHitList;
        pCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitList);

        CaloHitVector caloHitVector(caloHitList.begin(), caloHitList.end());
        std::sort(caloHitVector.begin(), caloHitVector.end(), LArClusterHelper::SortHitsByPosition);

        for (const CaloHit *const pCaloHit : caloHitVector)
        {
            if (pCaloHit->GetPositionVector().GetDotProduct(localWeightedDirection) < evtProjectedVtxPos)
            {
                minBeforeProjectedPos = std::min(minBeforeProjectedPos, pCaloHit->GetPositionVector().GetDotProduct(localWeightedDirection));
                maxBeforeProjectedPos = std::max(maxBeforeProjectedPos, pCaloHit->GetPositionVector().GetDotProduct(localWeightedDirection));

                beforeVtxEnergy += pCaloHit->GetElectromagneticEnergy();
                ++beforeVtxHits;
            }

            else
            {
                minAfterProjectedPos = std::min(minAfterProjectedPos, pCaloHit->GetPositionVector().GetDotProduct(localWeightedDirection));
                maxAfterProjectedPos = std::max(maxAfterProjectedPos, pCaloHit->GetPositionVector().GetDotProduct(localWeightedDirection));

                afterVtxEnergy += pCaloHit->GetElectromagneticEnergy();
                ++afterVtxHits;
            }
        }
    }

    // Use energy metrics if possible, otherwise fall back on hit counting.
    const unsigned int totalHits(beforeVtxHits + afterVtxHits);

    const float beforeLength(std::fabs(maxBeforeProjectedPos - minBeforeProjectedPos));
    const float afterLength(std::fabs(maxAfterProjectedPos - minAfterProjectedPos));

    const float beforeVtxEnergyDeposition(beforeLength < std::numeric_limits<float>::epsilon() ? 0 : beforeVtxEnergy / beforeLength);

    const float afterVtxEnergyDeposition(afterLength < std::numeric_limits<float>::epsilon() ? 0 : afterVtxEnergy / afterLength);

    const float totalEnergyDeposition(beforeVtxEnergyDeposition + afterVtxEnergyDeposition);

    if (useEnergyMetrics && totalEnergyDeposition > std::numeric_limits<float>::epsilon())
        return std::fabs((afterVtxEnergyDeposition - beforeVtxEnergyDeposition)) / totalEnergyDeposition;

    if (0 == totalHits)
        throw StatusCodeException(STATUS_CODE_FAILURE);

    const float beforeVtxHitDeposition(beforeLength < std::numeric_limits<float>::epsilon() ? 0 : beforeVtxHits / beforeLength);

    const float afterVtxHitDeposition(afterLength < std::numeric_limits<float>::epsilon() ? 0 : afterVtxHits / afterLength);

    const float totalHitDeposition(beforeVtxHitDeposition + afterVtxHitDeposition);

    if (totalHitDeposition > std::numeric_limits<float>::epsilon())
        return std::fabs((afterVtxHitDeposition - beforeVtxHitDeposition)) / totalHitDeposition;

    return 0.f;
}

//------------------------------------------------------------------------------------------------------------------------------------------

StatusCode EnergyDepositionAsymmetryFeatureTool::ReadSettings(const TiXmlHandle xmlHandle)
{
    return GlobalAsymmetryFeatureTool::ReadSettings(xmlHandle);
}

} // namespace lar_content