dd/da8/NeutrinoIdTool_8cc_source.html

 /**

  *  @file   larpandoracontent/LArControlFlow/NeutrinoIdTool.cc

  *

  *  @brief  Implementation of the neutrino id tool class.

  *

  *  $Log: $

  */


 #include "Pandora/AlgorithmHeaders.h"


 #include "larpandoracontent/LArControlFlow/NeutrinoIdTool.h"


 #include "Helpers/MCParticleHelper.h"


 #include "larpandoracontent/LArHelpers/LArFileHelper.h"

 #include "larpandoracontent/LArHelpers/LArMCParticleHelper.h"

 #include "larpandoracontent/LArHelpers/LArMvaHelper.h"

 #include "larpandoracontent/LArHelpers/LArPcaHelper.h"

 #include "larpandoracontent/LArHelpers/LArPfoHelper.h"


 #include "larpandoracontent/LArObjects/LArThreeDSlidingFitResult.h"


 using namespace pandora;


 namespace lar_content

 {


 template <typename T>

 NeutrinoIdTool<T>::NeutrinoIdTool() :

     m_useTrainingMode(false),

     m_selectNuanceCode(false),

     m_nuance(-std::numeric_limits<int>::max()),

     m_minPurity(0.9f),

     m_minCompleteness(0.9f),

     m_minProbability(0.0f),

     m_maxNeutrinos(1),

     m_persistFeatures(false),

     m_filePathEnvironmentVariable("FW_SEARCH_PATH")

 {

 }


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


 template <typename T>

 void NeutrinoIdTool<T>::SelectOutputPfos(const Algorithm *const pAlgorithm, const SliceHypotheses &nuSliceHypotheses,

     const SliceHypotheses &crSliceHypotheses, PfoList &selectedPfos)

 {

     if (nuSliceHypotheses.size() != crSliceHypotheses.size())

         throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);


     const unsigned int nSlices(nuSliceHypotheses.size());

     if (nSlices == 0)

         return;


     SliceFeaturesVector sliceFeaturesVector;

     this->GetSliceFeatures(this, nuSliceHypotheses, crSliceHypotheses, sliceFeaturesVector);


     if (m_useTrainingMode)

     {

         // ATTN in training mode, just return everything as a cosmic-ray

         this->SelectAllPfos(pAlgorithm, crSliceHypotheses, selectedPfos);


         unsigned int bestSliceIndex(std::numeric_limits<unsigned int>::max());

         if (!this->GetBestMCSliceIndex(pAlgorithm, nuSliceHypotheses, crSliceHypotheses, bestSliceIndex))

             return;


         for (unsigned int sliceIndex = 0; sliceIndex < nSlices; ++sliceIndex)

         {

             const SliceFeatures &features(sliceFeaturesVector.at(sliceIndex));

             if (!features.IsFeatureVectorAvailable())

                 continue;


             LArMvaHelper::MvaFeatureVector featureVector;

             features.GetFeatureVector(featureVector);

             LArMvaHelper::ProduceTrainingExample(m_trainingOutputFile, sliceIndex == bestSliceIndex, featureVector);

         }


         return;

     }


     this->SelectPfosByProbability(pAlgorithm, nuSliceHypotheses, crSliceHypotheses, sliceFeaturesVector, selectedPfos);

 }


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


 template <typename T>

 void NeutrinoIdTool<T>::GetSliceFeatures(const NeutrinoIdTool<T> *const pTool, const SliceHypotheses &nuSliceHypotheses,

     const SliceHypotheses &crSliceHypotheses, SliceFeaturesVector &sliceFeaturesVector) const

 {

     for (unsigned int sliceIndex = 0, nSlices = nuSliceHypotheses.size(); sliceIndex < nSlices; ++sliceIndex)

         sliceFeaturesVector.push_back(SliceFeatures(nuSliceHypotheses.at(sliceIndex), crSliceHypotheses.at(sliceIndex), pTool));

 }


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


 template <typename T>

 bool NeutrinoIdTool<T>::GetBestMCSliceIndex(const Algorithm *const pAlgorithm, const SliceHypotheses &nuSliceHypotheses,

     const SliceHypotheses &crSliceHypotheses, unsigned int &bestSliceIndex) const

 {

     unsigned int nHitsInBestSlice(0), nNuHitsInBestSlice(0);


     // Get all hits in all slices to find true number of mc hits

     const CaloHitList *pAllReconstructedCaloHitList(nullptr);

     PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*pAlgorithm, pAllReconstructedCaloHitList));


     const MCParticleList *pMCParticleList(nullptr);

     PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*pAlgorithm, pMCParticleList));


     // Obtain map: [mc particle -> primary mc particle]

     LArMCParticleHelper::MCRelationMap mcToPrimaryMCMap;

     LArMCParticleHelper::GetMCPrimaryMap(pMCParticleList, mcToPrimaryMCMap);


     // Remove non-reconstructable hits, e.g. those downstream of a neutron

     CaloHitList reconstructableCaloHitList;

     LArMCParticleHelper::PrimaryParameters parameters;

     LArMCParticleHelper::SelectCaloHits(pAllReconstructedCaloHitList, mcToPrimaryMCMap, reconstructableCaloHitList,

         parameters.m_selectInputHits, parameters.m_maxPhotonPropagation);


     const int nuNHitsTotal(this->CountNeutrinoInducedHits(reconstructableCaloHitList));

     const CaloHitSet reconstructableCaloHitSet(reconstructableCaloHitList.begin(), reconstructableCaloHitList.end());


     for (unsigned int sliceIndex = 0, nSlices = nuSliceHypotheses.size(); sliceIndex < nSlices; ++sliceIndex)

     {

         CaloHitList reconstructedCaloHitList;

         this->Collect2DHits(crSliceHypotheses.at(sliceIndex), reconstructedCaloHitList, reconstructableCaloHitSet);


         for (const ParticleFlowObject *const pNeutrino : nuSliceHypotheses.at(sliceIndex))

         {

             const PfoList &nuFinalStates(pNeutrino->GetDaughterPfoList());

             this->Collect2DHits(nuFinalStates, reconstructedCaloHitList, reconstructableCaloHitSet);

         }


         const unsigned int nNuHits(this->CountNeutrinoInducedHits(reconstructedCaloHitList));


         if (nNuHits > nNuHitsInBestSlice)

         {

             nNuHitsInBestSlice = nNuHits;

             nHitsInBestSlice = reconstructedCaloHitList.size();

             bestSliceIndex = sliceIndex;

         }

     }


     // ATTN for events with no neutrino induced hits, default neutrino purity and completeness to zero

     const float purity(nHitsInBestSlice > 0 ? static_cast<float>(nNuHitsInBestSlice) / static_cast<float>(nHitsInBestSlice) : 0.f);

     const float completeness(nuNHitsTotal > 0 ? static_cast<float>(nNuHitsInBestSlice) / static_cast<float>(nuNHitsTotal) : 0.f);

     return this->PassesQualityCuts(pAlgorithm, purity, completeness);

 }


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


 template <typename T>

 bool NeutrinoIdTool<T>::PassesQualityCuts(const Algorithm *const pAlgorithm, const float purity, const float completeness) const

 {

     if (purity < m_minPurity || completeness < m_minCompleteness)

         return false;

     if (m_selectNuanceCode && (this->GetNuanceCode(pAlgorithm) != m_nuance))

         return false;


     return true;

 }


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


 template <typename T>

 void NeutrinoIdTool<T>::Collect2DHits(const PfoList &pfos, CaloHitList &reconstructedCaloHitList, const CaloHitSet &reconstructableCaloHitSet) const

 {

     CaloHitList collectedHits;

     LArPfoHelper::GetCaloHits(pfos, TPC_VIEW_U, collectedHits);

     LArPfoHelper::GetCaloHits(pfos, TPC_VIEW_V, collectedHits);

     LArPfoHelper::GetCaloHits(pfos, TPC_VIEW_W, collectedHits);


     for (const CaloHit *const pCaloHit : collectedHits)

     {

         const CaloHit *const pParentHit = static_cast<const CaloHit *>(pCaloHit->GetParentAddress());

         if (!reconstructableCaloHitSet.count(pParentHit))

             continue;


         // Ensure no hits have been double counted

         if (std::find(reconstructedCaloHitList.begin(), reconstructedCaloHitList.end(), pParentHit) == reconstructedCaloHitList.end())

             reconstructedCaloHitList.push_back(pParentHit);

     }

 }


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


 template <typename T>

 unsigned int NeutrinoIdTool<T>::CountNeutrinoInducedHits(const CaloHitList &caloHitList) const

 {

     unsigned int nNuHits(0);

     for (const CaloHit *const pCaloHit : caloHitList)

     {

         try

         {

             if (LArMCParticleHelper::IsNeutrino(LArMCParticleHelper::GetParentMCParticle(MCParticleHelper::GetMainMCParticle(pCaloHit))))

                 nNuHits++;

         }

         catch (const StatusCodeException &)

         {

         }

     }


     return nNuHits;

 }


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


 template <typename T>

 int NeutrinoIdTool<T>::GetNuanceCode(const Algorithm *const pAlgorithm) const

 {

     const MCParticleList *pMCParticleList = nullptr;

     PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*pAlgorithm, pMCParticleList));


     MCParticleVector trueNeutrinos;

     LArMCParticleHelper::GetTrueNeutrinos(pMCParticleList, trueNeutrinos);


     if (trueNeutrinos.size() != 1)

     {

         std::cout << "NeutrinoIdTool::GetNuanceCode - Error: number of true neutrinos in event must be exactly one" << std::endl;

         throw StatusCodeException(STATUS_CODE_OUT_OF_RANGE);

     }


     return LArMCParticleHelper::GetNuanceCode(trueNeutrinos.front());

 }


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


 template <typename T>

 void NeutrinoIdTool<T>::SelectAllPfos(const pandora::Algorithm *const pAlgorithm, const SliceHypotheses &hypotheses, PfoList &selectedPfos) const

 {

     for (const PfoList &pfos : hypotheses)

     {

         for (const ParticleFlowObject *const pPfo : pfos)

         {

             object_creation::ParticleFlowObject::Metadata metadata;

             metadata.m_propertiesToAdd["NuScore"] = -1.f;

             PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::AlterMetadata(*pAlgorithm, pPfo, metadata));

         }


         this->SelectPfos(pfos, selectedPfos);

     }

 }


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


 template <typename T>

 void NeutrinoIdTool<T>::SelectPfosByProbability(const pandora::Algorithm *const pAlgorithm, const SliceHypotheses &nuSliceHypotheses,

     const SliceHypotheses &crSliceHypotheses, const SliceFeaturesVector &sliceFeaturesVector, PfoList &selectedPfos) const

 {

     // Calculate the probability of each slice that passes the minimum probability cut

     std::vector<UintFloatPair> sliceIndexProbabilityPairs;

     for (unsigned int sliceIndex = 0, nSlices = nuSliceHypotheses.size(); sliceIndex < nSlices; ++sliceIndex)

     {

         const float nuProbability(sliceFeaturesVector.at(sliceIndex).GetNeutrinoProbability(m_mva));


         for (const ParticleFlowObject *const pPfo : crSliceHypotheses.at(sliceIndex))

         {

             object_creation::ParticleFlowObject::Metadata metadata;

             metadata.m_propertiesToAdd["NuScore"] = nuProbability;


             if (m_persistFeatures)

             {

                 LArMvaHelper::DoubleMap featureMap;

                 sliceFeaturesVector.at(sliceIndex).GetFeatureMap(featureMap);


                 for (auto const &[name, value] : featureMap)

                     metadata.m_propertiesToAdd[name] = value;

             }


             PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::AlterMetadata(*pAlgorithm, pPfo, metadata));

         }


         for (const ParticleFlowObject *const pPfo : nuSliceHypotheses.at(sliceIndex))

         {

             object_creation::ParticleFlowObject::Metadata metadata;

             metadata.m_propertiesToAdd["NuScore"] = nuProbability;


             if (m_persistFeatures)

             {

                 LArMvaHelper::DoubleMap featureMap;

                 sliceFeaturesVector.at(sliceIndex).GetFeatureMap(featureMap);


                 for (auto const &[name, value] : featureMap)

                     metadata.m_propertiesToAdd[name] = value;

             }


             PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::AlterMetadata(*pAlgorithm, pPfo, metadata));

         }


         if (nuProbability < m_minProbability)

         {

             this->SelectPfos(crSliceHypotheses.at(sliceIndex), selectedPfos);

             continue;

         }


         sliceIndexProbabilityPairs.push_back(UintFloatPair(sliceIndex, nuProbability));

     }


     // Sort the slices by probability

     std::sort(sliceIndexProbabilityPairs.begin(), sliceIndexProbabilityPairs.end(),

         [](const UintFloatPair &a, const UintFloatPair &b) { return (a.second > b.second); });


     // Select the first m_maxNeutrinos as neutrinos, and the rest as cosmic

     unsigned int nNuSlices(0);

     for (const UintFloatPair &slice : sliceIndexProbabilityPairs)

     {

         if (nNuSlices < m_maxNeutrinos)

         {

             this->SelectPfos(nuSliceHypotheses.at(slice.first), selectedPfos);

             nNuSlices++;

             continue;

         }


         this->SelectPfos(crSliceHypotheses.at(slice.first), selectedPfos);

     }

 }


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


 template <typename T>

 void NeutrinoIdTool<T>::SelectPfos(const PfoList &pfos, PfoList &selectedPfos) const

 {

     selectedPfos.insert(selectedPfos.end(), pfos.begin(), pfos.end());

 }


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

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


 // TODO think about how to make this function cleaner when features are more established

 template <typename T>

 NeutrinoIdTool<T>::SliceFeatures::SliceFeatures(const PfoList &nuPfos, const PfoList &crPfos, const NeutrinoIdTool<T> *const pTool) :

     m_isAvailable(false),

     m_pTool(pTool)

 {

     try

     {

         const ParticleFlowObject *const pNeutrino(this->GetNeutrino(nuPfos));

         const CartesianVector &nuVertex(LArPfoHelper::GetVertex(pNeutrino)->GetPosition());

         const PfoList &nuFinalStates(pNeutrino->GetDaughterPfoList());


         // Neutrino features

         CartesianVector nuWeightedDirTotal(0.f, 0.f, 0.f);

         unsigned int nuNHitsUsedTotal(0);

         unsigned int nuNHitsTotal(0);

         CartesianPointVector nuAllSpacePoints;

         for (const ParticleFlowObject *const pPfo : nuFinalStates)

         {

             CartesianPointVector spacePoints;

             this->GetSpacePoints(pPfo, spacePoints);


             nuAllSpacePoints.insert(nuAllSpacePoints.end(), spacePoints.begin(), spacePoints.end());

             nuNHitsTotal += spacePoints.size();


             if (spacePoints.size() < 5)

                 continue;


             const CartesianVector dir(this->GetDirectionFromVertex(spacePoints, nuVertex));

             nuWeightedDirTotal += dir * static_cast<float>(spacePoints.size());

             nuNHitsUsedTotal += spacePoints.size();

         }


         if (nuNHitsUsedTotal == 0)

             return;

         const CartesianVector nuWeightedDir(nuWeightedDirTotal * (1.f / static_cast<float>(nuNHitsUsedTotal)));


         CartesianPointVector pointsInSphere;

         this->GetPointsInSphere(nuAllSpacePoints, nuVertex, 10, pointsInSphere);


         CartesianVector centroid(std::numeric_limits<float>::max(), std::numeric_limits<float>::max(), std::numeric_limits<float>::max());

         LArPcaHelper::EigenValues eigenValues(std::numeric_limits<float>::max(), std::numeric_limits<float>::max(), std::numeric_limits<float>::max());

         LArPcaHelper::EigenVectors eigenVectors;

         LArPcaHelper::RunPca(pointsInSphere, centroid, eigenValues, eigenVectors);


         const float nuNFinalStatePfos(static_cast<float>(nuFinalStates.size()));

         const float nuVertexY(nuVertex.GetY());

         const float nuWeightedDirZ(nuWeightedDir.GetZ());

         const float nuNSpacePointsInSphere(static_cast<float>(pointsInSphere.size()));


         if (eigenValues.GetX() <= std::numeric_limits<float>::epsilon())

             return;

         const float nuEigenRatioInSphere(eigenValues.GetY() / eigenValues.GetX());


         // Cosmic-ray features

         unsigned int nCRHitsMax(0);

         unsigned int nCRHitsTotal(0);

         float crLongestTrackDirY(std::numeric_limits<float>::max());

         float crLongestTrackDeflection(-std::numeric_limits<float>::max());


         for (const ParticleFlowObject *const pPfo : crPfos)

         {

             CartesianPointVector spacePoints;

             this->GetSpacePoints(pPfo, spacePoints);


             nCRHitsTotal += spacePoints.size();


             if (spacePoints.size() < 5)

                 continue;


             if (spacePoints.size() > nCRHitsMax)

             {

                 nCRHitsMax = spacePoints.size();

                 const CartesianVector upperDir(this->GetUpperDirection(spacePoints));

                 const CartesianVector lowerDir(this->GetLowerDirection(spacePoints));


                 crLongestTrackDirY = upperDir.GetY();

                 crLongestTrackDeflection = 1.f - upperDir.GetDotProduct(lowerDir * (-1.f));

             }

         }


         if (nCRHitsMax == 0)

             return;

         if (nCRHitsTotal == 0)

             return;


         const float crFracHitsInLongestTrack = static_cast<float>(nCRHitsMax) / static_cast<float>(nCRHitsTotal);


         // Push the features to the feature vector

         m_featureVector.push_back(nuNFinalStatePfos);

         m_featureVector.push_back(nuNHitsTotal);

         m_featureVector.push_back(nuVertexY);

         m_featureVector.push_back(nuWeightedDirZ);

         m_featureVector.push_back(nuNSpacePointsInSphere);

         m_featureVector.push_back(nuEigenRatioInSphere);

         m_featureVector.push_back(crLongestTrackDirY);

         m_featureVector.push_back(crLongestTrackDeflection);

         m_featureVector.push_back(crFracHitsInLongestTrack);

         m_featureVector.push_back(nCRHitsMax);


         m_featureMap["NuNFinalStatePfos"] = nuNFinalStatePfos;

         m_featureMap["NuNHitsTotal"] = nuNHitsTotal;

         m_featureMap["NuVertexY"] = nuVertexY;

         m_featureMap["NuWeightedDirZ"] = nuWeightedDirZ;

         m_featureMap["NuNSpacePointsInSphere"] = nuNSpacePointsInSphere;

         m_featureMap["NuEigenRatioInSphere"] = nuEigenRatioInSphere;

         m_featureMap["CRLongestTrackDirY"] = crLongestTrackDirY;

         m_featureMap["CRLongestTrackDeflection"] = crLongestTrackDeflection;

         m_featureMap["CRFracHitsInLongestTrack"] = crFracHitsInLongestTrack;

         m_featureMap["CRNHitsMax"] = nCRHitsMax;


         m_isAvailable = true;

     }

     catch (StatusCodeException &)

     {

         return;

     }

 }


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


 template <typename T>

 bool NeutrinoIdTool<T>::SliceFeatures::IsFeatureVectorAvailable() const

 {

     return m_isAvailable;

 }


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


 template <typename T>

 void NeutrinoIdTool<T>::SliceFeatures::GetFeatureVector(LArMvaHelper::MvaFeatureVector &featureVector) const

 {

     if (!m_isAvailable)

         throw StatusCodeException(STATUS_CODE_NOT_FOUND);


     featureVector.insert(featureVector.end(), m_featureVector.begin(), m_featureVector.end());

 }


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


 template <typename T>

 void NeutrinoIdTool<T>::SliceFeatures::GetFeatureMap(LArMvaHelper::DoubleMap &featureMap) const

 {

     if (!m_isAvailable)

         throw StatusCodeException(STATUS_CODE_NOT_FOUND);


     featureMap.insert(m_featureMap.begin(), m_featureMap.end());

 }


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


 template <typename T>

 float NeutrinoIdTool<T>::SliceFeatures::GetNeutrinoProbability(const T &t) const

 {

     // ATTN if one or more of the features can not be calculated, then default to calling the slice a cosmic ray

     if (!this->IsFeatureVectorAvailable())

         return 0.f;


     LArMvaHelper::MvaFeatureVector featureVector;

     this->GetFeatureVector(featureVector);

     return LArMvaHelper::CalculateProbability(t, featureVector);

 }


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


 template <typename T>

 const ParticleFlowObject *NeutrinoIdTool<T>::SliceFeatures::GetNeutrino(const PfoList &nuPfos) const

 {

     // ATTN we should only ever have one neutrino reconstructed per slice

     if (nuPfos.size() != 1)

         throw StatusCodeException(STATUS_CODE_OUT_OF_RANGE);


     return nuPfos.front();

 }


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


 template <typename T>

 void NeutrinoIdTool<T>::SliceFeatures::GetSpacePoints(const ParticleFlowObject *const pPfo, CartesianPointVector &spacePoints) const

 {

     ClusterList clusters3D;

     LArPfoHelper::GetThreeDClusterList(pPfo, clusters3D);


     if (clusters3D.size() > 1)

         throw StatusCodeException(STATUS_CODE_OUT_OF_RANGE);


     if (clusters3D.empty())

         return;


     CaloHitList caloHits;

     clusters3D.front()->GetOrderedCaloHitList().FillCaloHitList(caloHits);


     for (const CaloHit *const pCaloHit : caloHits)

         spacePoints.push_back(pCaloHit->GetPositionVector());

 }


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


 template <typename T>

 CartesianVector NeutrinoIdTool<T>::SliceFeatures::GetDirectionFromVertex(const CartesianPointVector &spacePoints, const CartesianVector &vertex) const

 {

     return this->GetDirection(spacePoints, [&](const CartesianVector &pointA, const CartesianVector &pointB) {

         return ((pointA - vertex).GetMagnitude() < (pointB - vertex).GetMagnitude());

     });

 }


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


 template <typename T>

 CartesianVector NeutrinoIdTool<T>::SliceFeatures::GetUpperDirection(const CartesianPointVector &spacePoints) const

 {

     return this->GetDirection(

         spacePoints, [&](const CartesianVector &pointA, const CartesianVector &pointB) { return (pointA.GetY() > pointB.GetY()); });

 }


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


 template <typename T>

 CartesianVector NeutrinoIdTool<T>::SliceFeatures::GetLowerDirection(const CartesianPointVector &spacePoints) const

 {

     return this->GetDirection(

         spacePoints, [&](const CartesianVector &pointA, const CartesianVector &pointB) { return (pointA.GetY() < pointB.GetY()); });

 }


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


 template <typename T>

 CartesianVector NeutrinoIdTool<T>::SliceFeatures::GetDirection(const CartesianPointVector &spacePoints,

     std::function<bool(const CartesianVector &pointA, const CartesianVector &pointB)> fShouldChooseA) const

 {

     // ATTN If wire w pitches vary between TPCs, exception will be raised in initialisation of lar pseudolayer plugin

     const LArTPC *const pFirstLArTPC(m_pTool->GetPandora().GetGeometry()->GetLArTPCMap().begin()->second);

     const float layerPitch(pFirstLArTPC->GetWirePitchW());


     const ThreeDSlidingFitResult fit(&spacePoints, 5, layerPitch);

     const CartesianVector endMin(fit.GetGlobalMinLayerPosition());

     const CartesianVector endMax(fit.GetGlobalMaxLayerPosition());

     const CartesianVector dirMin(fit.GetGlobalMinLayerDirection());

     const CartesianVector dirMax(fit.GetGlobalMaxLayerDirection());


     const bool isMinStart(fShouldChooseA(endMin, endMax));

     const CartesianVector startPoint(isMinStart ? endMin : endMax);

     const CartesianVector endPoint(isMinStart ? endMax : endMin);

     const CartesianVector startDir(isMinStart ? dirMin : dirMax);


     const bool shouldFlip((endPoint - startPoint).GetUnitVector().GetDotProduct(startDir) < 0.f);

     return (shouldFlip ? startDir * (-1.f) : startDir);

 }


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


 template <typename T>

 void NeutrinoIdTool<T>::SliceFeatures::GetPointsInSphere(const CartesianPointVector &spacePoints, const CartesianVector &vertex,

     const float radius, CartesianPointVector &spacePointsInSphere) const

 {

     for (const CartesianVector &point : spacePoints)

     {

         if ((point - vertex).GetMagnitudeSquared() <= radius * radius)

             spacePointsInSphere.push_back(point);

     }

 }


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


 template <typename T>

 StatusCode NeutrinoIdTool<T>::ReadSettings(const TiXmlHandle xmlHandle)

 {

     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "UseTrainingMode", m_useTrainingMode));


     if (m_useTrainingMode)

     {

         PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "TrainingOutputFileName", m_trainingOutputFile));

     }


     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinimumPurity", m_minPurity));


     PANDORA_RETURN_RESULT_IF_AND_IF(

         STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinimumCompleteness", m_minCompleteness));


     PANDORA_RETURN_RESULT_IF_AND_IF(

         STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SelectNuanceCode", m_selectNuanceCode));


     if (m_selectNuanceCode)

     {

         PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "NuanceCode", m_nuance));

     }


     PANDORA_RETURN_RESULT_IF_AND_IF(

         STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinimumNeutrinoProbability", m_minProbability));


     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaximumNeutrinos", m_maxNeutrinos));


     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "PersistFeatures", m_persistFeatures));


     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,

         XmlHelper::ReadValue(xmlHandle, "FilePathEnvironmentVariable", m_filePathEnvironmentVariable));


     if (!m_useTrainingMode)

     {

         std::string mvaName;

         PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "MvaName", mvaName));


         std::string mvaFileName;

         PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "MvaFileName", mvaFileName));


         const std::string fullMvaFileName(LArFileHelper::FindFileInPath(mvaFileName, m_filePathEnvironmentVariable));

         m_mva.Initialize(fullMvaFileName, mvaName);

     }


     return STATUS_CODE_SUCCESS;

 }


 template class NeutrinoIdTool<AdaBoostDecisionTree>;

 template class NeutrinoIdTool<SupportVectorMachine>;


 } // namespace lar_content

vertex
process_name vertex
Definition: cheaterreco.fcl:51

lar_content::NeutrinoIdTool::PassesQualityCuts
bool PassesQualityCuts(const pandora::Algorithm *const pAlgorithm, const float purity, const float completeness) const
Determine if the event passes the selection cuts for training and has the required NUANCE code...
Definition: NeutrinoIdTool.cc:152

pmtana::std
double std(const std::vector< short > &wf, const double ped_mean, size_t start, size_t nsample)
Definition: UtilFunc.cxx:42

LArPfoHelper.h
Header file for the pfo helper class.

lar_content::NeutrinoIdTool::m_mva
T m_mva
The mva.
Definition: NeutrinoIdTool.h:272

NeutrinoIdTool.h
Header file for the neutrino id tool class.

lar_content::LArMCParticleHelper::PrimaryParameters::m_selectInputHits
bool m_selectInputHits
whether to select input hits
Definition: LArMCParticleHelper.h:66

lar_content::LArMvaHelper::MvaFeatureVector
MvaTypes::MvaFeatureVector MvaFeatureVector
Definition: LArMvaHelper.h:72

lar_content::NeutrinoIdTool::GetBestMCSliceIndex
bool GetBestMCSliceIndex(const pandora::Algorithm *const pAlgorithm, const SliceHypotheses &nuSliceHypotheses, const SliceHypotheses &crSliceHypotheses, unsigned int &bestSliceIndex) const
Get the slice with the most neutrino induced hits using Monte-Carlo information.
Definition: NeutrinoIdTool.cc:97

lar_content::NeutrinoIdTool::SliceFeatures::SliceFeatures
SliceFeatures(const pandora::PfoList &nuPfos, const pandora::PfoList &crPfos, const NeutrinoIdTool *const pTool)
Constructor.
Definition: NeutrinoIdTool.cc:330

lar_content::LArMvaHelper::ProduceTrainingExample
static pandora::StatusCode ProduceTrainingExample(const std::string &trainingOutputFile, const bool result, TCONTAINER &&featureContainer)
Produce a training example with the given features and result.
Definition: LArMvaHelper.h:282

lar_content::NeutrinoIdTool::SliceFeatures::GetFeatureMap
void GetFeatureMap(LArMvaHelper::DoubleMap &featureMap) const
Get the feature map for the MVA.
Definition: NeutrinoIdTool.cc:469

lar_content::NeutrinoIdTool::SliceFeatures::IsFeatureVectorAvailable
bool IsFeatureVectorAvailable() const
Check if all features were calculable.
Definition: NeutrinoIdTool.cc:450

lar_content::NeutrinoIdTool::m_minProbability
float m_minProbability
Minimum probability required to classify a slice as the neutrino.
Definition: NeutrinoIdTool.h:267

lar_content::NeutrinoIdTool::GetNuanceCode
int GetNuanceCode(const pandora::Algorithm *const pAlgorithm) const
Use the current MCParticle list to get the nuance code of the neutrino in the event.
Definition: NeutrinoIdTool.cc:208

lar_content::LArPcaHelper::EigenValues
pandora::CartesianVector EigenValues
Definition: LArPcaHelper.h:24

lar_content::NeutrinoIdTool::SliceFeatures
Slice features class.
Definition: NeutrinoIdTool.h:47

lar_content::NeutrinoIdTool::ReadSettings
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
Definition: NeutrinoIdTool.cc:593

lar_content::LArPfoHelper::GetVertex
static const pandora::Vertex * GetVertex(const pandora::ParticleFlowObject *const pPfo)
Get the pfo vertex.
Definition: LArPfoHelper.cc:471

lar_content::NeutrinoIdTool::SelectOutputPfos
void SelectOutputPfos(const pandora::Algorithm *const pAlgorithm, const SliceHypotheses &nuSliceHypotheses, const SliceHypotheses &crSliceHypotheses, pandora::PfoList &selectedPfos)
Select which reconstruction hypotheses to use; neutrino outcomes or cosmic-ray muon outcomes for each...
Definition: NeutrinoIdTool.cc:45

lar_content::NeutrinoIdTool::SelectAllPfos
void SelectAllPfos(const pandora::Algorithm *const pAlgorithm, const SliceHypotheses &hypotheses, pandora::PfoList &selectedPfos) const
Select all pfos under the same hypothesis.
Definition: NeutrinoIdTool.cc:228

lar_content::NeutrinoIdTool::SliceFeatures::m_featureMap
LArMvaHelper::DoubleMap m_featureMap
A map between MVA features and their names.
Definition: NeutrinoIdTool.h:157

LArPcaHelper.h
Header file for the principal curve analysis helper class.

lar_content::LArMCParticleHelper::GetMCPrimaryMap
static void GetMCPrimaryMap(const pandora::MCParticleList *const pMCParticleList, MCRelationMap &mcPrimaryMap)
Get mapping from individual mc particles (in a provided list) and their primary parent mc particles...
Definition: LArMCParticleHelper.cc:418

lar_content::LArMCParticleHelper::GetNuanceCode
static unsigned int GetNuanceCode(const pandora::MCParticle *const pMCParticle)
Get the nuance code of an MCParticle.
Definition: LArMCParticleHelper.cc:119

lar_content::NeutrinoIdTool::m_minCompleteness
float m_minCompleteness
Minimum completeness of the best slice to use event for training.
Definition: NeutrinoIdTool.h:264

larg4::false
return false
Definition: OpFastScintillation.hh:592

lar_content::LArMCParticleHelper::PrimaryParameters::m_maxPhotonPropagation
float m_maxPhotonPropagation
the maximum photon propagation length
Definition: LArMCParticleHelper.h:67

a
process_name gaushit a
Definition: decode_signalprocess_icarus.fcl:37

lar_content::NeutrinoIdTool::m_filePathEnvironmentVariable
std::string m_filePathEnvironmentVariable
The environment variable providing a list of paths to mva files.
Definition: NeutrinoIdTool.h:273

lar_content::NeutrinoIdTool::Collect2DHits
void Collect2DHits(const pandora::PfoList &pfos, pandora::CaloHitList &reconstructedCaloHitList, const pandora::CaloHitSet &reconstructableCaloHitSet) const
Collect all 2D hits in a supplied list of Pfos and push them on to an existing hit list...
Definition: NeutrinoIdTool.cc:165

lar_content::NeutrinoIdTool::m_selectNuanceCode
bool m_selectNuanceCode
Should select training events by nuance code.
Definition: NeutrinoIdTool.h:261

lar_content::NeutrinoIdTool::m_maxNeutrinos
unsigned int m_maxNeutrinos
The maximum number of neutrinos to select in any one event.
Definition: NeutrinoIdTool.h:268

LArMvaHelper.h

lar_content::NeutrinoIdTool::CountNeutrinoInducedHits
unsigned int CountNeutrinoInducedHits(const pandora::CaloHitList &caloHitList) const
Count the number of neutrino induced hits in a given list using MC information.
Definition: NeutrinoIdTool.cc:187

lar_content::NeutrinoIdTool::UintFloatPair
std::pair< unsigned int, float > UintFloatPair
Definition: NeutrinoIdTool.h:161

lar_content::NeutrinoIdTool::SliceFeatures::GetFeatureVector
void GetFeatureVector(LArMvaHelper::MvaFeatureVector &featureVector) const
Get the feature vector for the MVA.
Definition: NeutrinoIdTool.cc:458

lar_content::NeutrinoIdTool::SliceFeatures::GetNeutrinoProbability
float GetNeutrinoProbability(const T &t) const
Get the probability that this slice contains a neutrino interaction.
Definition: NeutrinoIdTool.cc:480

lar_content::NeutrinoIdTool::SliceFeatures::GetLowerDirection
pandora::CartesianVector GetLowerDirection(const pandora::CartesianPointVector &spacePoints) const
Use a sliding fit to get the lower direction of a collection of spacepoints.
Definition: NeutrinoIdTool.cc:546

lar_content::NeutrinoIdTool::GetSliceFeatures
void GetSliceFeatures(const NeutrinoIdTool *const pTool, const SliceHypotheses &nuSliceHypotheses, const SliceHypotheses &crSliceHypotheses, SliceFeaturesVector &sliceFeaturesVector) const
Get the features of each slice.
Definition: NeutrinoIdTool.cc:87

lar_content::LArMvaHelper::DoubleMap
std::map< std::string, double > DoubleMap
Definition: LArMvaHelper.h:73

LArMCParticleHelper.h
Header file for the lar monte carlo particle helper helper class.

lar_content::NeutrinoIdTool::m_nuance
int m_nuance
Nuance code to select for training.
Definition: NeutrinoIdTool.h:262

lar_content::NeutrinoIdTool::SliceFeatures::m_isAvailable
bool m_isAvailable
Is the feature vector available.
Definition: NeutrinoIdTool.h:155

lar_content::NeutrinoIdTool::SelectPfos
void SelectPfos(const pandora::PfoList &pfos, pandora::PfoList &selectedPfos) const
Add the given pfos to the selected Pfo list.
Definition: NeutrinoIdTool.cc:320

lar_content::NeutrinoIdTool::m_useTrainingMode
bool m_useTrainingMode
Should use training mode. If true, training examples will be written to the output file...
Definition: NeutrinoIdTool.h:259

lar_content::SliceHypotheses
std::vector< pandora::PfoList > SliceHypotheses
Definition: MasterAlgorithm.h:29

lar_content::NeutrinoIdTool::m_minPurity
float m_minPurity
Minimum purity of the best slice to use event for training.
Definition: NeutrinoIdTool.h:263

lar_content::NeutrinoIdTool::SliceFeatures::m_featureVector
LArMvaHelper::MvaFeatureVector m_featureVector
The MVA feature vector.
Definition: NeutrinoIdTool.h:156

lar_content::NeutrinoIdTool::SliceFeatures::GetNeutrino
const pandora::ParticleFlowObject * GetNeutrino(const pandora::PfoList &nuPfos) const
Get the recontructed neutrino the input list of neutrino Pfos.
Definition: NeutrinoIdTool.cc:494

lar_pandora::MCParticleVector
std::vector< art::Ptr< simb::MCParticle > > MCParticleVector
Definition: LArPandoraHelper.h:67

lar_content::LArMCParticleHelper::GetParentMCParticle
static const pandora::MCParticle * GetParentMCParticle(const pandora::MCParticle *const pMCParticle)
Get the parent mc particle.
Definition: LArMCParticleHelper.cc:232

lar_content::LArFileHelper::FindFileInPath
static std::string FindFileInPath(const std::string &unqualifiedFileName, const std::string &environmentVariable, const std::string &delimiter=":")
Find the fully-qualified file name by searching through a list of delimiter-separated paths in a name...
Definition: LArFileHelper.cc:21

LArFileHelper.h
Header file for the file helper class.

lar_content::LArPfoHelper::GetThreeDClusterList
static void GetThreeDClusterList(const pandora::ParticleFlowObject *const pPfo, pandora::ClusterList &clusterList)
Get the list of 3D clusters from an input pfo.
Definition: LArPfoHelper.cc:136

lar_content::NeutrinoIdTool::SelectPfosByProbability
void SelectPfosByProbability(const pandora::Algorithm *const pAlgorithm, const SliceHypotheses &nuSliceHypotheses, const SliceHypotheses &crSliceHypotheses, const SliceFeaturesVector &sliceFeaturesVector, pandora::PfoList &selectedPfos) const
Select pfos based on the probability that their slice contains a neutrino interaction.
Definition: NeutrinoIdTool.cc:246

lar_content::NeutrinoIdTool::SliceFeaturesVector
std::vector< SliceFeatures > SliceFeaturesVector
Definition: NeutrinoIdTool.h:162

LArThreeDSlidingFitResult.h
Header file for the lar three dimensional sliding fit result class.

dropbox.dir
tuple dir
Definition: dropbox.py:28

lar_content::NeutrinoIdTool::m_persistFeatures
bool m_persistFeatures
If true, the mva features will be persisted in the metadata.
Definition: NeutrinoIdTool.h:270

lar_content::NeutrinoIdTool
NeutrinoIdTool class.
Definition: NeutrinoIdTool.h:32

lar_content::NeutrinoIdTool::SliceFeatures::GetDirectionFromVertex
pandora::CartesianVector GetDirectionFromVertex(const pandora::CartesianPointVector &spacePoints, const pandora::CartesianVector &vertex) const
Use a sliding fit to get the direction of a collection of spacepoint near a vertex position...
Definition: NeutrinoIdTool.cc:527

lar_content::NeutrinoIdTool::SliceFeatures::GetUpperDirection
pandora::CartesianVector GetUpperDirection(const pandora::CartesianPointVector &spacePoints) const
Use a sliding fit to get the upper direction of a collection of spacepoints.
Definition: NeutrinoIdTool.cc:537

lar_content::NeutrinoIdTool::SliceFeatures::GetDirection
pandora::CartesianVector GetDirection(const pandora::CartesianPointVector &spacePoints, std::function< bool(const pandora::CartesianVector &pointA, const pandora::CartesianVector &pointB)> fShouldChooseA) const
Use a sliding fit to get the direction of a collection of spacepoints.
Definition: NeutrinoIdTool.cc:555

lar_content::LArPcaHelper::RunPca
static void RunPca(const T &t, pandora::CartesianVector &centroid, EigenValues &outputEigenValues, EigenVectors &outputEigenVectors)
Run principal component analysis using input calo hits (TPC_VIEW_U,V,W or TPC_3D; all treated as 3D p...

pandora
required by fuzzyCluster table::sbnd_g4_services gaushitTruthMatch pandora
Definition: reco_sbnd.fcl:182

lar_content::LArMCParticleHelper::GetTrueNeutrinos
static void GetTrueNeutrinos(const pandora::MCParticleList *const pMCParticleList, pandora::MCParticleVector &trueNeutrinos)
Get neutrino MC particles from an input MC particle list.
Definition: LArMCParticleHelper.cc:206

lar_content::LArPcaHelper::EigenVectors
std::vector< pandora::CartesianVector > EigenVectors
Definition: LArPcaHelper.h:25

lar_content::ThreeDSlidingFitResult
ThreeDSlidingFitResult class.
Definition: LArThreeDSlidingFitResult.h:23

name
then echo fcl name
Definition: sbndpoms_genfclwithrunnumber_maker.sh:220

value
temporary value
Definition: pmtsimulation_icarus.fcl:216

lar_content::NeutrinoIdTool::SliceFeatures::GetPointsInSphere
void GetPointsInSphere(const pandora::CartesianPointVector &spacePoints, const pandora::CartesianVector &vertex, const float radius, pandora::CartesianPointVector &spacePointsInSphere) const
Get a vector of spacepoints within a given radius of a vertex point.
Definition: NeutrinoIdTool.cc:580

lar_content::NeutrinoIdTool::SliceFeatures::GetSpacePoints
void GetSpacePoints(const pandora::ParticleFlowObject *const pPfo, pandora::CartesianPointVector &spacePoints) const
Get the 3D space points in a given pfo.
Definition: NeutrinoIdTool.cc:506

lar_content::NeutrinoIdTool::m_trainingOutputFile
std::string m_trainingOutputFile
Output file name for training examples.
Definition: NeutrinoIdTool.h:260

lar_content::LArMvaHelper::CalculateProbability
static double CalculateProbability(const MvaInterface &classifier, TCONTAINER &&featureContainer)
Use the trained mva to calculate a classification probability for an example.
Definition: LArMvaHelper.h:366

lar_content::LArPfoHelper::GetCaloHits
static void GetCaloHits(const pandora::PfoList &pfoList, const pandora::HitType &hitType, pandora::CaloHitList &caloHitList)
Get a list of calo hits of a particular hit type from a list of pfos.

lar_content::LArMCParticleHelper::MCRelationMap
std::unordered_map< const pandora::MCParticle *, const pandora::MCParticle * > MCRelationMap
Definition: LArMCParticleHelper.h:27

cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10

lar_content::LArMCParticleHelper::IsNeutrino
static bool IsNeutrino(const pandora::MCParticle *const pMCParticle)
Whether a mc particle is a neutrino or antineutrino.
Definition: LArMCParticleHelper.cc:131

lar_content::LArMCParticleHelper::PrimaryParameters
PrimaryParameters class.
Definition: LArMCParticleHelper.h:55

lar_content::LArMCParticleHelper::SelectCaloHits
static void SelectCaloHits(const pandora::CaloHitList *const pCaloHitList, const MCRelationMap &mcToTargetMCMap, pandora::CaloHitList &selectedCaloHitList, const bool selectInputHits, const float maxPhotonPropagation)
Select a subset of calo hits representing those that represent &quot;reconstructable&quot; regions of the event...
Definition: LArMCParticleHelper.cc:1116