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