LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
VertexMonitoringAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
12 
18 
19 using namespace pandora;
20 
21 namespace lar_content
22 {
23 
24 VertexMonitoringAlgorithm::VertexMonitoringAlgorithm() :
25  m_visualise{true},
26  m_writeFile{false},
29  m_scalingFactor{1.f}
30 {
31 }
32 
33 //------------------------------------------------------------------------------------------------------------------------------------------
34 
36 {
37  if (m_writeFile)
38  {
39  PANDORA_MONITORING_API(SaveTree(this->GetPandora(), m_treename.c_str(), m_filename.c_str(), "UPDATE"));
40  }
41 }
42 
43 //------------------------------------------------------------------------------------------------------------------------------------------
44 
46 {
47  if (m_visualise)
48  {
49  PANDORA_MONITORING_API(SetEveDisplayParameters(
50  this->GetPandora(), true, DETECTOR_VIEW_XZ, m_transparencyThresholdE, m_energyScaleThresholdE, m_scalingFactor));
51  }
52 
53  this->AssessVertices();
54 
55  if (m_visualise)
56  {
57  PANDORA_MONITORING_API(ViewEvent(this->GetPandora()));
58  }
59 
60  return STATUS_CODE_SUCCESS;
61 }
62 
63 //------------------------------------------------------------------------------------------------------------------------------------------
64 
66 {
67 #ifdef MONITORING
68  const MCParticleList *pMCParticleList{nullptr};
69  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pMCParticleList));
70  const PfoList *pPfoList{nullptr};
71  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pPfoList));
72 
74  MCParticleVector primaries;
75  LArMCParticleHelper::GetPrimaryMCParticleList(pMCParticleList, primaries);
76  const MCParticle *pTrueNeutrino{nullptr};
77  const ParticleFlowObject *pRecoNeutrino{nullptr};
78  if (!primaries.empty())
79  {
80  for (const MCParticle *primary : primaries)
81  {
82  const MCParticleList &parents{primary->GetParentList()};
83  if (parents.size() == 1 && LArMCParticleHelper::IsNeutrino(parents.front()))
84  {
85  pTrueNeutrino = parents.front();
86  break;
87  }
88  }
89  }
90 
91  for (const ParticleFlowObject *pPfo : *pPfoList)
92  {
93  if (LArPfoHelper::IsNeutrino(pPfo))
94  {
95  pRecoNeutrino = pPfo;
96  break;
97  }
98  }
99 
100  MCParticleList primariesList(primaries.begin(), primaries.end());
102 
103  if (pRecoNeutrino && pTrueNeutrino)
104  {
105  const LArTransformationPlugin *transform{this->GetPandora().GetPlugins()->GetLArTransformationPlugin()};
106  const CartesianVector &trueVertex{pTrueNeutrino->GetVertex()};
107  const CartesianVector &recoVertex{LArPfoHelper::GetVertex(pRecoNeutrino)->GetPosition()};
108  if (m_visualise)
109  {
110  const CartesianVector tu(trueVertex.GetX(), 0.f, static_cast<float>(transform->YZtoU(trueVertex.GetY(), trueVertex.GetZ())));
111  const CartesianVector tv(trueVertex.GetX(), 0.f, static_cast<float>(transform->YZtoV(trueVertex.GetY(), trueVertex.GetZ())));
112  const CartesianVector tw(trueVertex.GetX(), 0.f, static_cast<float>(transform->YZtoW(trueVertex.GetY(), trueVertex.GetZ())));
113 
114  const CartesianVector ru(recoVertex.GetX(), 0.f, static_cast<float>(transform->YZtoU(recoVertex.GetY(), recoVertex.GetZ())));
115  const CartesianVector rv(recoVertex.GetX(), 0.f, static_cast<float>(transform->YZtoV(recoVertex.GetY(), recoVertex.GetZ())));
116  const CartesianVector rw(recoVertex.GetX(), 0.f, static_cast<float>(transform->YZtoW(recoVertex.GetY(), recoVertex.GetZ())));
117 
118  const float du{(ru - tu).GetMagnitude()};
119  const float dv{(rv - tv).GetMagnitude()};
120  const float dw{(rw - tw).GetMagnitude()};
121 
122  std::cout << "delta(u, v, w): (" << du << ", " << dv << "," << dw << ")" << std::endl;
123 
124  PANDORA_MONITORING_API(AddMarkerToVisualization(this->GetPandora(), &tu, "U true vertex", BLUE, 2));
125  PANDORA_MONITORING_API(AddMarkerToVisualization(this->GetPandora(), &tv, "V true vertex", BLUE, 2));
126  PANDORA_MONITORING_API(AddMarkerToVisualization(this->GetPandora(), &tw, "W true vertex", BLUE, 2));
127  PANDORA_MONITORING_API(AddMarkerToVisualization(this->GetPandora(), &ru, "U reco vertex", RED, 2));
128  PANDORA_MONITORING_API(AddMarkerToVisualization(this->GetPandora(), &rv, "V reco vertex", RED, 2));
129  PANDORA_MONITORING_API(AddMarkerToVisualization(this->GetPandora(), &rw, "W reco vertex", RED, 2));
130  }
131 
132  if (m_writeFile && LArVertexHelper::IsInFiducialVolume(this->GetPandora(), trueVertex, "dune_fd_hd"))
133  {
134  const CartesianVector delta{recoVertex - trueVertex};
135  const float dx{delta.GetX()}, dy{delta.GetY()}, dz{delta.GetZ()}, dr{delta.GetMagnitude()};
136  const float trueNuEnergy{pTrueNeutrino->GetEnergy()};
137  const int success{1};
138  const int isCC{descriptor.IsCC()};
139  const int isQE{descriptor.IsQE()};
140  const int isRes{descriptor.IsResonant()};
141  const int isDIS{descriptor.IsDIS()};
142  const int isCoh{descriptor.IsCoherent()};
143  const int isOther{!(isCC || isQE || isRes || isDIS)};
144  const int isMu{descriptor.IsMuonNeutrino()};
145  const int isElectron{descriptor.IsElectronNeutrino()};
146  const int nPiZero{static_cast<int>(descriptor.GetNumPiZero())};
147  const int nPiPlus{static_cast<int>(descriptor.GetNumPiPlus())};
148  const int nPiMinus{static_cast<int>(descriptor.GetNumPiMinus())};
149  const int nPhotons{static_cast<int>(descriptor.GetNumPhotons())};
150  const int nProtons{static_cast<int>(descriptor.GetNumProtons())};
151  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "success", success));
152  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "trueNuEnergy", trueNuEnergy));
153  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isCC", isCC));
154  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isQE", isQE));
155  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isRes", isRes));
156  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isDIS", isDIS));
157  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isCoh", isCoh));
158  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isOther", isOther));
159  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isMu", isMu));
160  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isElectron", isElectron));
161  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "nPiZero", nPiZero));
162  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "nPiPlus", nPiPlus));
163  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "nPiMinus", nPiMinus));
164  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "nPhotons", nPhotons));
165  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "nProtons", nProtons));
166  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "dx", dx));
167  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "dy", dy));
168  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "dz", dz));
169  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "dr", dr));
170  PANDORA_MONITORING_API(FillTree(this->GetPandora(), m_treename.c_str()));
171  }
172  }
173  else if (pTrueNeutrino)
174  {
175  const CartesianVector &trueVertex{pTrueNeutrino->GetVertex()};
176 
177  if (m_writeFile && LArVertexHelper::IsInFiducialVolume(this->GetPandora(), trueVertex, "dune_fd_hd"))
178  {
179  const int success{0};
180  const float dx{-999.f}, dy{-999.f}, dz{-999.f}, dr{-999.f};
181  const float trueNuEnergy{pTrueNeutrino->GetEnergy()};
182  const int isCC{descriptor.IsCC()};
183  const int isQE{descriptor.IsQE()};
184  const int isRes{descriptor.IsResonant()};
185  const int isDIS{descriptor.IsDIS()};
186  const int isCoh{descriptor.IsCoherent()};
187  const int isOther{!(isCC || isQE || isRes || isDIS)};
188  const int isMu{descriptor.IsMuonNeutrino()};
189  const int isElectron{descriptor.IsElectronNeutrino()};
190  const int nPiZero{static_cast<int>(descriptor.GetNumPiZero())};
191  const int nPiPlus{static_cast<int>(descriptor.GetNumPiPlus())};
192  const int nPiMinus{static_cast<int>(descriptor.GetNumPiMinus())};
193  const int nPhotons{static_cast<int>(descriptor.GetNumPhotons())};
194  const int nProtons{static_cast<int>(descriptor.GetNumProtons())};
195  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "success", success));
196  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "trueNuEnergy", trueNuEnergy));
197  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isCC", isCC));
198  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isQE", isQE));
199  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isRes", isRes));
200  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isDIS", isDIS));
201  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isCoh", isCoh));
202  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isOther", isOther));
203  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isMu", isMu));
204  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "isElectron", isElectron));
205  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "nPiZero", nPiZero));
206  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "nPiPlus", nPiPlus));
207  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "nPiMinus", nPiMinus));
208  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "nPhotons", nPhotons));
209  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "nProtons", nProtons));
210  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "dx", dx));
211  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "dy", dy));
212  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "dz", dz));
213  PANDORA_MONITORING_API(SetTreeVariable(this->GetPandora(), m_treename.c_str(), "dr", dr));
214  PANDORA_MONITORING_API(FillTree(this->GetPandora(), m_treename.c_str()));
215  }
216  }
217 #endif
218 
219  return STATUS_CODE_SUCCESS;
220 }
221 
222 //------------------------------------------------------------------------------------------------------------------------------------------
223 
224 StatusCode VertexMonitoringAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
225 {
226  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "Visualize", m_visualise));
227  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "WriteFile", m_writeFile));
228  if (m_writeFile)
229  {
230  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "Filename", m_filename));
231  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "Treename", m_treename));
232  }
233 
234  PANDORA_RETURN_RESULT_IF_AND_IF(
235  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "TransparencyThresholdE", m_transparencyThresholdE));
236  PANDORA_RETURN_RESULT_IF_AND_IF(
237  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "EnergyScaleThresholdE", m_energyScaleThresholdE));
238  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ScalingFactor", m_scalingFactor));
239 
240  return STATUS_CODE_SUCCESS;
241 }
242 
243 } // namespace lar_content
Header file for the pfo helper class.
std::unordered_map< const pandora::MCParticle *, pandora::CaloHitList > MCContributionMap
Header file for the interaction type helper class.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
static const pandora::Vertex * GetVertex(const pandora::ParticleFlowObject *const pPfo)
Get the pfo vertex.
float m_energyScaleThresholdE
Cell energy for which color is at top end of continous color palette.
float m_scalingFactor
TEve works with [cm], Pandora usually works with [mm] (but LArContent went with cm too) ...
Header file for the geometry helper class.
Header file for the lar monte carlo particle helper helper class.
std::vector< art::Ptr< simb::MCParticle > > MCParticleVector
static bool IsNeutrino(const pandora::ParticleFlowObject *const pPfo)
Whether a pfo is a neutrino or (antineutrino)
static bool IsInFiducialVolume(const pandora::Pandora &pandora, const pandora::CartesianVector &vertex, const std::string &detector)
Determine if a vertex is within a detector&#39;s fiducial volume. This throws a STATUS_CODE_INVALID_PARAM...
Header file for the vertex helper class.
static void GetPrimaryMCParticleList(const pandora::MCParticleList *const pMCParticleList, pandora::MCParticleVector &mcPrimaryVector)
Get vector of primary MC particles from an input list of MC particles.
Header file for the particle visualisation algorithm.
static InteractionDescriptor GetInteractionDescriptor(const pandora::MCParticleList &mcPrimaryList)
Get the interaction descriptor of an event.
float m_transparencyThresholdE
Cell energy for which transparency is saturated (0%, fully opaque)
static bool IsNeutrino(const pandora::MCParticle *const pMCParticle)
Whether a mc particle is a neutrino or antineutrino.