LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
CheatingNeutrinoCreationAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
13 
15 
16 using namespace pandora;
17 
18 namespace lar_content
19 {
20 
21 CheatingNeutrinoCreationAlgorithm::CheatingNeutrinoCreationAlgorithm() :
22  m_collapseToPrimaryMCParticles(false),
23  m_vertexTolerance(0.5f)
24 {
25 }
26 
27 //------------------------------------------------------------------------------------------------------------------------------------------
28 
30 {
31  MCParticleVector mcNeutrinoVector;
32  this->GetMCNeutrinoVector(mcNeutrinoVector);
33 
34  for (const MCParticle *const pMCNeutrino : mcNeutrinoVector)
35  {
36  const ParticleFlowObject *pNeutrinoPfo(nullptr);
37  this->CreateAndSaveNeutrinoPfo(pMCNeutrino, pNeutrinoPfo);
38 
39  if (!pNeutrinoPfo)
40  return STATUS_CODE_FAILURE;
41 
42  this->AddNeutrinoVertex(pMCNeutrino, pNeutrinoPfo);
43 
44  MCParticleToPfoMap mcParticleToPfoMap;
45  this->GetMCParticleToDaughterPfoMap(mcParticleToPfoMap);
46  this->CreatePfoHierarchy(pMCNeutrino, pNeutrinoPfo, mcParticleToPfoMap);
47  }
48 
49  return STATUS_CODE_SUCCESS;
50 }
51 
52 //------------------------------------------------------------------------------------------------------------------------------------------
53 
55 {
56  const MCParticleList *pMCParticleList(nullptr);
57  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetList(*this, m_mcParticleListName, pMCParticleList));
58 
59  MCParticleVector allMCNeutrinoVector;
60  LArMCParticleHelper::GetTrueNeutrinos(pMCParticleList, allMCNeutrinoVector);
61 
62  for (const MCParticle *const pMCNeutrino : allMCNeutrinoVector)
63  {
64  if (LArMCParticleHelper::IsNeutrino(pMCNeutrino) && (MC_3D == pMCNeutrino->GetMCParticleType()))
65  mcNeutrinoVector.push_back(pMCNeutrino);
66  }
67 }
68 
69 //------------------------------------------------------------------------------------------------------------------------------------------
70 
71 void CheatingNeutrinoCreationAlgorithm::CreateAndSaveNeutrinoPfo(const MCParticle *const pMCNeutrino, const ParticleFlowObject *&pNeutrinoPfo) const
72 {
73  pNeutrinoPfo = nullptr;
74 
75  PandoraContentApi::ParticleFlowObject::Parameters pfoParameters;
76  pfoParameters.m_particleId = pMCNeutrino->GetParticleId();
77  pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
78  pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
79  pfoParameters.m_energy = pMCNeutrino->GetEnergy();
80  pfoParameters.m_momentum = pMCNeutrino->GetMomentum();
81 
82  std::string neutrinoPfoListName;
83  const PfoList *pNeutrinoPfoList(nullptr);
84  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*this, pNeutrinoPfoList, neutrinoPfoListName));
85  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*this, pfoParameters, pNeutrinoPfo));
86 
87  if (!pNeutrinoPfoList || pNeutrinoPfoList->empty())
88  throw StatusCodeException(STATUS_CODE_FAILURE);
89 
90  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Pfo>(*this, m_neutrinoPfoListName));
91  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ReplaceCurrentList<Pfo>(*this, m_neutrinoPfoListName));
92 }
93 
94 //------------------------------------------------------------------------------------------------------------------------------------------
95 
96 void CheatingNeutrinoCreationAlgorithm::AddNeutrinoVertex(const MCParticle *const pMCNeutrino, const ParticleFlowObject *const pNeutrinoPfo) const
97 {
98  const VertexList *pVertexList(nullptr);
99  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pVertexList));
100 
101  const Vertex *pNeutrinoVertex(nullptr);
102  float closestVertexDistance(std::numeric_limits<float>::max());
103 
104  for (const Vertex *const pVertex : *pVertexList)
105  {
106  const float distance((pVertex->GetPosition() - pMCNeutrino->GetEndpoint()).GetMagnitude());
107 
108  if (distance < closestVertexDistance)
109  {
110  pNeutrinoVertex = pVertex;
111  closestVertexDistance = distance;
112  }
113  }
114 
115  if (!pNeutrinoVertex || (VERTEX_3D != pNeutrinoVertex->GetVertexType()) || ((pNeutrinoVertex->GetPosition() - pMCNeutrino->GetEndpoint()).GetMagnitude() > m_vertexTolerance))
116  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
117 
118  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToPfo(*this, pNeutrinoPfo, pNeutrinoVertex));
119 }
120 
121 //------------------------------------------------------------------------------------------------------------------------------------------
122 
124 {
126 
128  {
129  const MCParticleList *pMCParticleList(nullptr);
130  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetList(*this, m_mcParticleListName, pMCParticleList));
131 
132  LArMCParticleHelper::GetMCPrimaryMap(pMCParticleList, mcPrimaryMap);
133  }
134 
135  for (const std::string &daughterPfoListName : m_daughterPfoListNames)
136  {
137  const PfoList *pDaughterPfoList(nullptr);
138 
139  if (STATUS_CODE_SUCCESS != PandoraContentApi::GetList(*this, daughterPfoListName, pDaughterPfoList))
140  {
141  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
142  std::cout << "CheatingNeutrinoCreationAlgorithm: pfo list " << daughterPfoListName << " unavailable." << std::endl;
143 
144  continue;
145  }
146 
147  for (const ParticleFlowObject *const pDaughterPfo : *pDaughterPfoList)
148  {
149  const MCParticle *pMCParticle(LArMCParticleHelper::GetMainMCParticle(pDaughterPfo));
150 
152  {
153  LArMCParticleHelper::MCRelationMap::const_iterator primaryIter = mcPrimaryMap.find(pMCParticle);
154 
155  if (mcPrimaryMap.end() == primaryIter)
156  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
157 
158  pMCParticle = primaryIter->second;
159  }
160 
161  if (!mcParticleToPfoMap.insert(MCParticleToPfoMap::value_type(pMCParticle, pDaughterPfo)).second)
162  throw StatusCodeException(STATUS_CODE_OUT_OF_RANGE);
163  }
164  }
165 }
166 
167 //------------------------------------------------------------------------------------------------------------------------------------------
168 
169 void CheatingNeutrinoCreationAlgorithm::CreatePfoHierarchy(const MCParticle *const pParentMCParticle, const ParticleFlowObject *const pParentPfo,
170  const MCParticleToPfoMap &mcParticleToPfoMap) const
171 {
172  for (const MCParticle *const pDaughterMCParticle : pParentMCParticle->GetDaughterList())
173  {
174  MCParticleToPfoMap::const_iterator mapIter = mcParticleToPfoMap.find(pDaughterMCParticle);
175 
176  if (mcParticleToPfoMap.end() == mapIter)
177  {
178  this->CreatePfoHierarchy(pDaughterMCParticle, pParentPfo, mcParticleToPfoMap);
179  continue;
180  }
181 
182  const ParticleFlowObject *const pDaughterPfo(mapIter->second);
183  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SetPfoParentDaughterRelationship(*this, pParentPfo, pDaughterPfo));
184  this->CreatePfoHierarchy(pDaughterMCParticle, pDaughterPfo, mcParticleToPfoMap);
185  }
186 }
187 
188 //------------------------------------------------------------------------------------------------------------------------------------------
189 
190 StatusCode CheatingNeutrinoCreationAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
191 {
192  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
193  "CollapseToPrimaryMCParticles", m_collapseToPrimaryMCParticles));
194 
195  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
196  "MCParticleListName", m_mcParticleListName));
197 
198  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
199  "NeutrinoPfoListName", m_neutrinoPfoListName));
200 
201  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
202  "VertexListName", m_vertexListName));
203 
204  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadVectorOfValues(xmlHandle,
205  "DaughterPfoListNames", m_daughterPfoListNames));
206 
207  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
208  "VertexTolerance", m_vertexTolerance));
209 
210  return STATUS_CODE_SUCCESS;
211 }
212 
213 } // namespace lar_content
Header file for the pfo helper class.
std::string m_neutrinoPfoListName
The name of the neutrino pfo list.
float m_vertexTolerance
Tolerance, 3d displacement, allowed between reco neutrino vertex and mc neutrino endpoint.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
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...
static const pandora::MCParticle * GetMainMCParticle(const pandora::ParticleFlowObject *const pPfo)
Find the mc particle making the largest contribution to 2D clusters in a specified pfo...
std::vector< art::Ptr< simb::MCParticle > > MCParticleVector
std::unordered_map< const pandora::MCParticle *, const pandora::ParticleFlowObject * > MCParticleToPfoMap
TFile f
Definition: plotHisto.C:6
std::string m_mcParticleListName
The name of the three d mc particle list name.
void CreatePfoHierarchy(const pandora::MCParticle *const pParentMCParticle, const pandora::ParticleFlowObject *const pParentPfo, const MCParticleToPfoMap &mcParticleToPfoMap) const
Use information from mc particles and the mc particle to pfo map to fully-reconstruct the daughter pf...
Int_t max
Definition: plot.C:27
intermediate_table::const_iterator const_iterator
Header file for the lar monte carlo particle helper helper class.
pandora::StringVector m_daughterPfoListNames
The list of daughter pfo list names.
void GetMCParticleToDaughterPfoMap(MCParticleToPfoMap &mcParticleToPfoMap) const
Extract candidate daughter pfos from external lists and populate a map from main mc particle (or prim...
static void GetTrueNeutrinos(const pandora::MCParticleList *const pMCParticleList, pandora::MCParticleVector &trueNeutrinos)
Get neutrino MC particles from an input MC particle list.
Header file for the cheating neutrino creation algorithm class.
std::string m_vertexListName
The name of the neutrino vertex list.
bool m_collapseToPrimaryMCParticles
Whether to collapse mc particle hierarchies to primary particles.
void GetMCNeutrinoVector(pandora::MCParticleVector &mcNeutrinoVector) const
Get the mc neutrino vector.
std::unordered_map< const pandora::MCParticle *, const pandora::MCParticle * > MCRelationMap
void AddNeutrinoVertex(const pandora::MCParticle *const pMCNeutrino, const pandora::ParticleFlowObject *const pNeutrinoPfo) const
Extract reconstructed vertex from external list, check its position agrees with mc neutrino...
std::list< Vertex > VertexList
Definition: DCEL.h:178
static bool IsNeutrino(const pandora::MCParticle *const pMCParticle)
Whether a mc particle is a neutrino or antineutrino.
void CreateAndSaveNeutrinoPfo(const pandora::MCParticle *const pMCNeutrino, const pandora::ParticleFlowObject *&pNeutrinoPfo) const
Create and save a neutrino pfo with properties dictated by the mc neutrino.