ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RootPlanarClusterWriter.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RootPlanarClusterWriter.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include "ACTFW/Io/Root/RootPlanarClusterWriter.hpp"
10 
11 #include <TFile.h>
12 #include <TTree.h>
13 #include <ios>
14 #include <stdexcept>
15 
16 #include "ACTFW/EventData/SimHit.hpp"
17 #include "ACTFW/EventData/SimIdentifier.hpp"
18 #include "ACTFW/EventData/SimParticle.hpp"
19 #include "ACTFW/Framework/WhiteBoard.hpp"
20 #include "ACTFW/Utilities/Paths.hpp"
21 #include "Acts/Plugins/Digitization/DigitizationModule.hpp"
22 #include "Acts/Plugins/Digitization/PlanarModuleCluster.hpp"
23 #include "Acts/Plugins/Digitization/Segmentation.hpp"
24 #include "Acts/Plugins/Identification/IdentifiedDetectorElement.hpp"
25 #include "Acts/Utilities/Units.hpp"
26 
27 FW::RootPlanarClusterWriter::RootPlanarClusterWriter(
28  const FW::RootPlanarClusterWriter::Config& cfg, Acts::Logging::Level lvl)
29  : WriterT(cfg.inputClusters, "RootPlanarClusterWriter", lvl),
30  m_cfg(cfg),
31  m_outputFile(cfg.rootFile) {
32  // inputClusters is already checked by base constructor
33  if (m_cfg.inputSimulatedHits.empty()) {
34  throw std::invalid_argument("Missing simulated hits input collection");
35  }
36  if (m_cfg.treeName.empty()) {
37  throw std::invalid_argument("Missing tree name");
38  }
39  // Setup ROOT I/O
40  if (m_outputFile == nullptr) {
41  m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
42  if (m_outputFile == nullptr) {
43  throw std::ios_base::failure("Could not open '" + m_cfg.filePath);
44  }
45  }
46  m_outputFile->cd();
47  m_outputTree =
48  new TTree(m_cfg.treeName.c_str(), "TTree from RootPlanarClusterWriter");
49  if (m_outputTree == nullptr)
50  throw std::bad_alloc();
51 
52  // Set the branches
53  m_outputTree->Branch("event_nr", &m_eventNr);
54  m_outputTree->Branch("volume_id", &m_volumeID);
55  m_outputTree->Branch("layer_id", &m_layerID);
56  m_outputTree->Branch("surface_id", &m_surfaceID);
57  m_outputTree->Branch("g_x", &m_x);
58  m_outputTree->Branch("g_y", &m_y);
59  m_outputTree->Branch("g_z", &m_z);
60  m_outputTree->Branch("g_t", &m_t);
61  m_outputTree->Branch("l_x", &m_lx);
62  m_outputTree->Branch("l_y", &m_ly);
63  m_outputTree->Branch("cov_l_x", &m_cov_lx);
64  m_outputTree->Branch("cov_l_y", &m_cov_ly);
65  m_outputTree->Branch("cell_ID_x", &m_cell_IDx);
66  m_outputTree->Branch("cell_ID_y", &m_cell_IDy);
67  m_outputTree->Branch("cell_l_x", &m_cell_lx);
68  m_outputTree->Branch("cell_l_y", &m_cell_ly);
69  m_outputTree->Branch("cell_data", &m_cell_data);
70  m_outputTree->Branch("truth_g_x", &m_t_gx);
71  m_outputTree->Branch("truth_g_y", &m_t_gy);
72  m_outputTree->Branch("truth_g_z", &m_t_gz);
73  m_outputTree->Branch("truth_g_t", &m_t_gt);
74  m_outputTree->Branch("truth_l_x", &m_t_lx);
75  m_outputTree->Branch("truth_l_y", &m_t_ly);
76  m_outputTree->Branch("truth_barcode", &m_t_barcode, "truth_barcode/l");
77 }
78 
79 FW::RootPlanarClusterWriter::~RootPlanarClusterWriter() {
81  if (m_cfg.rootFile == nullptr) {
82  m_outputFile->Close();
83  }
84 }
85 
86 FW::ProcessCode FW::RootPlanarClusterWriter::endRun() {
87  // Write the tree
88  m_outputFile->cd();
89  m_outputTree->Write();
90  ACTS_INFO("Wrote particles to tree '" << m_cfg.treeName << "' in '"
91  << m_cfg.filePath << "'");
92  return ProcessCode::SUCCESS;
93 }
94 
95 FW::ProcessCode FW::RootPlanarClusterWriter::writeT(
96  const AlgorithmContext& ctx,
97  const FW::GeometryIdMultimap<Acts::PlanarModuleCluster>& clusters) {
98  // retrieve simulated hits
99  const auto& simHits =
100  ctx.eventStore.get<SimHitContainer>(m_cfg.inputSimulatedHits);
101 
102  // Exclusive access to the tree while writing
103  std::lock_guard<std::mutex> lock(m_writeMutex);
104  // Get the event number
105  m_eventNr = ctx.eventNumber;
106 
107  // Loop over the planar clusters in this event
108  for (const auto& entry : clusters) {
109  Acts::GeometryID geoId = entry.first;
110  const Acts::PlanarModuleCluster& cluster = entry.second;
111  // local cluster information: position, @todo coveraiance
112  auto parameters = cluster.parameters();
113  Acts::Vector2D local(parameters[Acts::ParDef::eLOC_0],
114  parameters[Acts::ParDef::eLOC_1]);
115 
117  Acts::Vector3D pos(0, 0, 0);
118  Acts::Vector3D mom(1, 1, 1);
119  // the cluster surface
120  const auto& clusterSurface = cluster.referenceSurface();
121  // transform local into global position information
122  clusterSurface.localToGlobal(ctx.geoContext, local, mom, pos);
123  // identification
124  m_volumeID = geoId.volume();
125  m_layerID = geoId.layer();
126  m_surfaceID = geoId.sensitive();
127  m_x = pos.x();
128  m_y = pos.y();
129  m_z = pos.z();
130  m_t = parameters[2] / Acts::UnitConstants::ns;
131  m_lx = local.x();
132  m_ly = local.y();
133  m_cov_lx = 0.; // @todo fill in
134  m_cov_ly = 0.; // @todo fill in
135  // get the cells and run through them
136  const auto& cells = cluster.digitizationCells();
137  auto detectorElement = dynamic_cast<const Acts::IdentifiedDetectorElement*>(
138  clusterSurface.associatedDetectorElement());
139  for (auto& cell : cells) {
140  // cell identification
141  m_cell_IDx.push_back(cell.channel0);
142  m_cell_IDy.push_back(cell.channel1);
143  m_cell_data.push_back(cell.data);
144  // for more we need the digitization module
145  if (detectorElement && detectorElement->digitizationModule()) {
146  auto digitationModule = detectorElement->digitizationModule();
147  const Acts::Segmentation& segmentation =
148  digitationModule->segmentation();
149  // get the cell positions
150  auto cellLocalPosition = segmentation.cellPosition(cell);
151  m_cell_lx.push_back(cellLocalPosition.x());
152  m_cell_ly.push_back(cellLocalPosition.y());
153  }
154  }
155  // write hit-particle truth association
156  // each hit can have multiple particles, e.g. in a dense environment
157  for (auto idx : cluster.sourceLink().indices()) {
158  auto it = simHits.nth(idx);
159  if (it == simHits.end()) {
160  ACTS_FATAL("Simulation hit with index " << idx << " does not exist");
161  return ProcessCode::ABORT;
162  }
163  const auto& simHit = *it;
164 
165  // local position to be calculated
166  Acts::Vector2D lPosition;
167  clusterSurface.globalToLocal(ctx.geoContext, simHit.position(),
168  simHit.unitDirection(), lPosition);
169  // fill the variables
170  m_t_gx.push_back(simHit.position().x());
171  m_t_gy.push_back(simHit.position().y());
172  m_t_gz.push_back(simHit.position().z());
173  m_t_gt.push_back(simHit.time());
174  m_t_lx.push_back(lPosition.x());
175  m_t_ly.push_back(lPosition.y());
176  m_t_barcode.push_back(simHit.particleId().value());
177  }
178  // fill the tree
179  m_outputTree->Fill();
180  // now reset
181  m_cell_IDx.clear();
182  m_cell_IDy.clear();
183  m_cell_lx.clear();
184  m_cell_ly.clear();
185  m_cell_data.clear();
186  m_t_gx.clear();
187  m_t_gy.clear();
188  m_t_gz.clear();
189  m_t_gt.clear();
190  m_t_lx.clear();
191  m_t_ly.clear();
192  m_t_barcode.clear();
193  }
194  return FW::ProcessCode::SUCCESS;
195 }