LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
EventWritingAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
11 #include "Persistency/BinaryFileWriter.h"
12 #include "Persistency/XmlFileWriter.h"
13 
16 
19 
21 
22 using namespace pandora;
23 
24 namespace lar_content
25 {
26 
27 EventWritingAlgorithm::EventWritingAlgorithm() :
28  m_geometryFileType(UNKNOWN_FILE_TYPE),
29  m_eventFileType(UNKNOWN_FILE_TYPE),
30  m_pEventFileWriter(nullptr),
31  m_pGeometryFileWriter(nullptr),
32  m_shouldWriteGeometry(false),
33  m_writtenGeometry(false),
34  m_shouldWriteEvents(true),
35  m_shouldWriteMCRelationships(true),
36  m_shouldWriteTrackRelationships(true),
37  m_shouldOverwriteEventFile(false),
38  m_shouldOverwriteGeometryFile(false),
39  m_useLArCaloHits(true),
40  m_useLArMCParticles(true),
41  m_shouldFilterByNuanceCode(false),
42  m_filterNuanceCode(0),
43  m_shouldFilterByMCParticles(false),
44  m_neutrinoInducedOnly(false),
45  m_matchingMinPrimaryHits(15),
46  m_nNonNeutrons(0),
47  m_nMuons(0),
48  m_nElectrons(0),
49  m_nProtons(0),
50  m_nPhotons(0),
51  m_nChargedPions(0),
52  m_shouldFilterByNeutrinoVertex(false),
53  m_detectorHalfLengthX(-1.f),
54  m_detectorHalfLengthY(-1.f),
55  m_detectorHalfLengthZ(-1.f),
56  m_coordinateOffsetX(0.f),
57  m_coordinateOffsetY(0.f),
58  m_coordinateOffsetZ(0.f),
59  m_selectedBorderX(-1.f),
60  m_selectedBorderY(-1.f),
61  m_selectedBorderZ(-1.f)
62 {
63 }
64 
65 //------------------------------------------------------------------------------------------------------------------------------------------
66 
68 {
69  delete m_pEventFileWriter;
70  delete m_pGeometryFileWriter;
71 }
72 
73 //------------------------------------------------------------------------------------------------------------------------------------------
74 
76 {
78  {
79  const FileMode fileMode(m_shouldOverwriteGeometryFile ? OVERWRITE : APPEND);
80 
81  if (BINARY == m_geometryFileType)
82  {
83  m_pGeometryFileWriter = new BinaryFileWriter(this->GetPandora(), m_geometryFileName, fileMode);
84  }
85  else if (XML == m_geometryFileType)
86  {
87  m_pGeometryFileWriter = new XmlFileWriter(this->GetPandora(), m_geometryFileName, fileMode);
88  }
89  else
90  {
91  return STATUS_CODE_FAILURE;
92  }
93  }
94 
96  {
97  const FileMode fileMode(m_shouldOverwriteEventFile ? OVERWRITE : APPEND);
98 
99  if (BINARY == m_eventFileType)
100  {
101  m_pEventFileWriter = new BinaryFileWriter(this->GetPandora(), m_eventFileName, fileMode);
102  }
103  else if (XML == m_eventFileType)
104  {
105  m_pEventFileWriter = new XmlFileWriter(this->GetPandora(), m_eventFileName, fileMode);
106  }
107  else
108  {
109  return STATUS_CODE_FAILURE;
110  }
111 
112  if (m_useLArCaloHits)
113  m_pEventFileWriter->SetFactory(new LArCaloHitFactory);
114 
116  m_pEventFileWriter->SetFactory(new LArMCParticleFactory);
117  }
118 
119  return STATUS_CODE_SUCCESS;
120 }
121 
122 //------------------------------------------------------------------------------------------------------------------------------------------
123 
125 {
126  // ATTN Should complete geometry creation in LArSoft begin job, but some channel status service functionality unavailable at that point
128  {
129  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, m_pGeometryFileWriter->WriteGeometry());
130  m_writtenGeometry = true;
131  }
132 
133  bool matchNuanceCode(!m_shouldFilterByNuanceCode || this->PassNuanceCodeFilter());
134  bool matchParticles(!m_shouldFilterByMCParticles || this->PassMCParticleFilter());
135  bool matchNeutrinoVertexPosition(!m_shouldFilterByNeutrinoVertex || this->PassNeutrinoVertexFilter());
136 
137  if (matchNuanceCode && matchParticles && matchNeutrinoVertexPosition && m_pEventFileWriter && m_shouldWriteEvents)
138  {
139  const CaloHitList *pCaloHitList = nullptr;
140  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pCaloHitList));
141 
142  const TrackList *pTrackList = nullptr;
143  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pTrackList));
144 
145  const MCParticleList *pMCParticleList = nullptr;
146  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pMCParticleList));
147 
148  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, m_pEventFileWriter->WriteEvent(*pCaloHitList, *pTrackList, *pMCParticleList,
150  }
151 
152  return STATUS_CODE_SUCCESS;
153 }
154 
155 //------------------------------------------------------------------------------------------------------------------------------------------
156 
158 {
159  const MCParticleList *pMCParticleList = nullptr;
160  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pMCParticleList));
161 
162  MCParticleVector mcNeutrinoList;
163  LArMCParticleHelper::GetTrueNeutrinos(pMCParticleList, mcNeutrinoList);
164 
165  for (const MCParticle *const pMCNeutrino : mcNeutrinoList)
166  {
167  const LArMCParticle *const pLArMCNeutrino = dynamic_cast<const LArMCParticle*>(pMCNeutrino);
168 
169  if (pLArMCNeutrino && (pLArMCNeutrino->GetNuanceCode() == m_filterNuanceCode))
170  return true;
171  }
172 
173  return false;
174 }
175 
176 //------------------------------------------------------------------------------------------------------------------------------------------
177 
179 {
180  const MCParticleList *pMCParticleList = nullptr;
181  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pMCParticleList));
182 
183  const CaloHitList *pCaloHitList = nullptr;
184  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pCaloHitList));
185 
187  LArMCParticleHelper::MCContributionMap mcParticlesToGoodHitsMap;
188  LArMCParticleHelper::SelectReconstructableMCParticles(pMCParticleList, pCaloHitList, parameters, LArMCParticleHelper::IsBeamNeutrinoFinalState, mcParticlesToGoodHitsMap);
189 
191  {
192  LArMCParticleHelper::SelectReconstructableMCParticles(pMCParticleList, pCaloHitList, parameters, LArMCParticleHelper::IsBeamParticle, mcParticlesToGoodHitsMap);
193  LArMCParticleHelper::SelectReconstructableMCParticles(pMCParticleList, pCaloHitList, parameters, LArMCParticleHelper::IsCosmicRay, mcParticlesToGoodHitsMap);
194  }
195 
196  unsigned int nNonNeutrons(0), nMuons(0), nElectrons(0), nProtons(0), nPhotons(0), nChargedPions(0);
197 
198  MCParticleVector mcPrimaryVector;
199  for (const auto &mapEntry : mcParticlesToGoodHitsMap) mcPrimaryVector.push_back(mapEntry.first);
200  std::sort(mcPrimaryVector.begin(), mcPrimaryVector.end(), LArMCParticleHelper::SortByMomentum);
201 
202  for (const MCParticle *const pMCPrimary : mcPrimaryVector)
203  {
204  const unsigned int particleId(std::abs(pMCPrimary->GetParticleId()));
205  if (NEUTRON != particleId) ++nNonNeutrons;
206 
207  if (MU_MINUS == particleId) ++nMuons;
208  else if (E_MINUS == particleId) ++nElectrons;
209  else if (PROTON == particleId) ++nProtons;
210  else if (PHOTON == particleId) ++nPhotons;
211  else if (PI_PLUS == particleId) ++nChargedPions;
212  }
213 
214  if ((nNonNeutrons == m_nNonNeutrons) && (nMuons == m_nMuons) && (nElectrons == m_nElectrons) &&
215  (nProtons == m_nProtons) && (nPhotons == m_nPhotons) && (nChargedPions == m_nChargedPions))
216  {
217  return true;
218  }
219 
220  return false;
221 }
222 
223 //------------------------------------------------------------------------------------------------------------------------------------------
224 
226 {
227  const MCParticleList *pMCParticleList = nullptr;
228  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pMCParticleList));
229 
230  MCParticleVector mcNeutrinoVector;
231  LArMCParticleHelper::GetTrueNeutrinos(pMCParticleList, mcNeutrinoVector);
232 
233  for (const MCParticle *const pMCNeutrino : mcNeutrinoVector)
234  {
235  const CartesianVector &neutrinoInteractionVertex(pMCNeutrino->GetEndpoint());
236 
237  if ((neutrinoInteractionVertex.GetX() < (m_detectorHalfLengthX - m_coordinateOffsetX - m_selectedBorderX)) && (neutrinoInteractionVertex.GetX() > (-m_coordinateOffsetX + m_selectedBorderX)) &&
238  (neutrinoInteractionVertex.GetY() < (m_detectorHalfLengthY - m_coordinateOffsetY - m_selectedBorderY)) && (neutrinoInteractionVertex.GetY() > (-m_coordinateOffsetY + m_selectedBorderY)) &&
239  (neutrinoInteractionVertex.GetZ() < (m_detectorHalfLengthZ - m_coordinateOffsetZ - m_selectedBorderZ)) && (neutrinoInteractionVertex.GetZ() > (-m_coordinateOffsetZ + m_selectedBorderZ)) )
240  {
241  return true;
242  }
243  }
244 
245  return false;
246 }
247 
248 //------------------------------------------------------------------------------------------------------------------------------------------
249 
250 StatusCode EventWritingAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
251 {
252  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
253  "ShouldWriteGeometry", m_shouldWriteGeometry));
254 
256  {
257  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
258  "GeometryFileName", m_geometryFileName));
259 
260  std::string fileExtension(m_geometryFileName.substr(m_geometryFileName.find_last_of(".")));
261  std::transform(fileExtension.begin(), fileExtension.end(), fileExtension.begin(), ::tolower);
262 
263  if (std::string(".xml") == fileExtension)
264  {
265  m_geometryFileType = XML;
266  }
267  else if (std::string(".pndr") == fileExtension)
268  {
269  m_geometryFileType = BINARY;
270  }
271  else
272  {
273  std::cout << "EventReadingAlgorithm: Unknown geometry file type specified " << std::endl;
274  return STATUS_CODE_INVALID_PARAMETER;
275  }
276  }
277 
278  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
279  "ShouldWriteEvents", m_shouldWriteEvents));
280 
282  {
283  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
284  "EventFileName", m_eventFileName));
285 
286  std::string fileExtension(m_eventFileName.substr(m_eventFileName.find_last_of(".")));
287  std::transform(fileExtension.begin(), fileExtension.end(), fileExtension.begin(), ::tolower);
288 
289  if (std::string(".xml") == fileExtension)
290  {
291  m_eventFileType = XML;
292  }
293  else if (std::string(".pndr") == fileExtension)
294  {
295  m_eventFileType = BINARY;
296  }
297  else
298  {
299  std::cout << "EventReadingAlgorithm: Unknown event file type specified " << std::endl;
300  return STATUS_CODE_INVALID_PARAMETER;
301  }
302  }
303 
304  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
305  "ShouldWriteMCRelationships", m_shouldWriteMCRelationships));
306 
307  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
308  "ShouldWriteTrackRelationships", m_shouldWriteTrackRelationships));
309 
310  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
311  "ShouldOverwriteEventFile", m_shouldOverwriteEventFile));
312 
313  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
314  "ShouldOverwriteGeometryFile", m_shouldOverwriteGeometryFile));
315 
316  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
317  "UseLArCaloHits", m_useLArCaloHits));
318 
319  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
320  "UseLArMCParticles", m_useLArMCParticles));
321 
322  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
323  "ShouldFilterByNuanceCode", m_shouldFilterByNuanceCode));
324 
325  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
326  "ShouldFilterByMCParticles", m_shouldFilterByMCParticles));
327 
328  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
329  "ShouldFilterByNeutrinoVertex", m_shouldFilterByNeutrinoVertex));
330 
332  {
333  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "FilterNuanceCode", m_filterNuanceCode));
334  }
335 
337  {
338  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "NeutrinoInducedOnly", m_neutrinoInducedOnly));
339  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "MatchingMinPrimaryHits", m_matchingMinPrimaryHits));
340  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "NNonNeutrons", m_nNonNeutrons));
341  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "NMuons", m_nMuons));
342  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "NElectrons", m_nElectrons));
343  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "NProtons", m_nProtons));
344  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "NPhotons", m_nPhotons));
345  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "NChargedPions", m_nChargedPions));
346  }
347 
349  {
350  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "DetectorHalfLengthX", m_detectorHalfLengthX));
351  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "DetectorHalfLengthY", m_detectorHalfLengthY));
352  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "DetectorHalfLengthZ", m_detectorHalfLengthZ));
353  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "CoordinateOffsetX", m_coordinateOffsetX));
354  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "CoordinateOffsetY", m_coordinateOffsetY));
355  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "CoordinateOffsetZ", m_coordinateOffsetZ));
356  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "SelectedBorderX", m_selectedBorderX));
357  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "SelectedBorderY", m_selectedBorderY));
358  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "SelectedBorderZ", m_selectedBorderZ));
359  }
360 
361  return STATUS_CODE_SUCCESS;
362 }
363 
364 } // namespace lar_content
float m_detectorHalfLengthZ
Half length of detector in z dimension.
unsigned int m_matchingMinPrimaryHits
The minimum number of mc primary hits used in matching scheme.
bool m_shouldWriteMCRelationships
Whether to write mc relationship information to the events file.
pandora::FileWriter * m_pEventFileWriter
Address of the event file writer.
Header file for the event writing algorithm class.
bool m_shouldFilterByMCParticles
Whether to filter output by mc particle constituents.
bool m_shouldWriteEvents
Whether to write events to a specified file.
Header file for the lar calo hit class.
bool PassNeutrinoVertexFilter() const
Whether current event passes neutrino vertex position filter (e.g. fiducial volume cut) ...
bool m_useLArCaloHits
Whether to write lar calo hits, or standard pandora calo hits.
float m_coordinateOffsetY
Origin offset (from detector corner) in y dimension.
int GetNuanceCode() const
Get the nuance code.
bool m_neutrinoInducedOnly
Whether to consider only mc particles that were neutrino induced.
unsigned int m_nPhotons
The requested number of mc primaries that are photons.
unsigned int m_nProtons
The requested number of mc primaries that are protons.
float m_selectedBorderY
Required distance from detector edge in y dimension.
float m_detectorHalfLengthY
Half length of detector in y dimension.
float m_selectedBorderZ
Required distance from detector edge in z dimension.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
std::string m_eventFileName
Name of the output event file.
Header file for the lar monitoring helper helper class.
std::vector< art::Ptr< simb::MCParticle > > MCParticleVector
pandora::FileType m_geometryFileType
The geometry file type.
unsigned int m_nChargedPions
The requested number of mc primaries that are charged pions.
TFile f
Definition: plotHisto.C:6
bool PassNuanceCodeFilter() const
Whether current event passes nuance code filter.
bool m_shouldWriteTrackRelationships
Whether to write track relationship information to the events file.
bool m_useLArMCParticles
Whether to write lar mc particles, or standard pandora mc particles.
static bool SortByMomentum(const pandora::MCParticle *const pLhs, const pandora::MCParticle *const pRhs)
Sort mc particles by their momentum.
LAr mc particle class.
Definition: LArMCParticle.h:38
static bool IsCosmicRay(const pandora::MCParticle *const pMCParticle)
Return true if passed a primary cosmic ray MCParticle.
pandora::FileWriter * m_pGeometryFileWriter
Address of the geometry file writer.
static void SelectReconstructableMCParticles(const pandora::MCParticleList *pMCParticleList, const pandora::CaloHitList *pCaloHitList, const PrimaryParameters &parameters, std::function< bool(const pandora::MCParticle *const)> fCriteria, MCContributionMap &selectedMCParticlesToHitsMap)
Select primary, reconstructable mc particles that match given criteria.
std::unordered_map< const pandora::MCParticle *, pandora::CaloHitList > MCContributionMap
unsigned int m_nNonNeutrons
The requested number of mc primaries that are not neutrons.
Header file for the lar monte carlo particle helper helper class.
static bool IsBeamParticle(const pandora::MCParticle *const pMCParticle)
Returns true if passed a primary beam MCParticle.
float m_selectedBorderX
Required distance from detector edge in x dimension.
bool m_shouldFilterByNeutrinoVertex
Whether to filter output by neutrino vertex position (e.g. fiducial volume cut)
bool PassMCParticleFilter() const
Whether current event passes mc particle constituent filter.
int m_filterNuanceCode
The filter nuance code (required if specify filter by nuance code)
float m_detectorHalfLengthX
Half length of detector in x dimension.
Header file for the lar mc particle class.
unsigned int m_nElectrons
The requested number of mc primaries that are electrons.
LArMCParticleFactory responsible for object creation.
Definition: LArMCParticle.h:64
std::string m_geometryFileName
Name of the output geometry file.
static void GetTrueNeutrinos(const pandora::MCParticleList *const pMCParticleList, pandora::MCParticleVector &trueNeutrinos)
Get neutrino MC particles from an input MC particle list.
pandora::FileType m_eventFileType
The event file type.
bool m_shouldFilterByNuanceCode
Whether to filter output by nuance code.
bool m_shouldOverwriteEventFile
Whether to overwrite existing event file with specified name, or append.
unsigned int m_nMuons
The requested number of mc primaries that are muons.
LArCaloHitFactory responsible for object creation.
Definition: LArCaloHit.h:64
bool m_shouldWriteGeometry
Whether to write geometry to a specified file.
bool m_shouldOverwriteGeometryFile
Whether to overwrite existing geometry file with specified name, or append.
float m_coordinateOffsetZ
Origin offset (from detector corner) in z dimension.
static bool IsBeamNeutrinoFinalState(const pandora::MCParticle *const pMCParticle)
Returns true if passed a primary neutrino final state MCParticle.
bool m_writtenGeometry
Whether geometry has been written.
float m_coordinateOffsetX
Origin offset (from detector corner) in x dimension.