Geant4HepMC/src/EventRecording.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Geant4HepMC/EventRecording.hpp"
0010
0011 #include "ActsExamples/EventData/SimParticle.hpp"
0012 #include "ActsExamples/Framework/WhiteBoard.hpp"
0013
0014 #include <stdexcept>
0015
0016 #include <FTFP_BERT.hh>
0017 #include <G4RunManager.hh>
0018 #include <G4VUserDetectorConstruction.hh>
0019 #include <HepMC3/GenParticle.h>
0020
0021 #include "EventAction.hpp"
0022 #include "PrimaryGeneratorAction.hpp"
0023 #include "RunAction.hpp"
0024 #include "SteppingAction.hpp"
0025
0026 namespace ActsExamples {
0027
0028 EventRecording::~EventRecording() {
0029   m_runManager = nullptr;
0030 }
0031
0032 EventRecording::EventRecording(const EventRecording::Config& config,
0033                                Acts::Logging::Level level)
0034     : IAlgorithm("EventRecording", level),
0035       m_cfg(config),
0036       m_runManager(std::make_unique<G4RunManager>()) {
0037   if (m_cfg.inputParticles.empty()) {
0038     throw std::invalid_argument("Missing input particle collection");
0039   }
0040   if (m_cfg.outputHepMcTracks.empty()) {
0041     throw std::invalid_argument("Missing output event collection");
0042   }
0043   if (m_cfg.detector == nullptr) {
0044     throw std::invalid_argument("Missing detector construction object");
0045   }
0046
0047   m_inputParticles.initialize(m_cfg.inputParticles);
0048   m_outputEvents.initialize(m_cfg.outputHepMcTracks);
0049
0050   // Now set up the Geant4 simulation
0051
0052   // G4RunManager deals with the lifetime of these objects
0053   m_runManager->SetUserInitialization(
0054       m_cfg.detector->buildGeant4DetectorConstruction(m_cfg.constructionOptions)
0055           .release());
0056   m_runManager->SetUserInitialization(new FTFP_BERT);
0057   m_runManager->SetUserAction(new Geant4::HepMC3::RunAction());
0058   m_runManager->SetUserAction(
0059       new Geant4::HepMC3::EventAction(m_cfg.processesCombine));
0060   m_runManager->SetUserAction(
0061       new Geant4::HepMC3::PrimaryGeneratorAction(m_cfg.seed1, m_cfg.seed2));
0062   m_runManager->SetUserAction(
0063       new Geant4::HepMC3::SteppingAction(m_cfg.processesReject));
0064   m_runManager->Initialize();
0065 }
0066
0067 ProcessCode EventRecording::execute(const AlgorithmContext& context) const {
0068   // ensure exclusive access to the geant run manager
0069   std::lock_guard<std::mutex> guard(m_runManagerLock);
0070
0071   // Retrieve the initial particles
0072   const auto initialParticles = m_inputParticles(context);
0073
0074   // Storage of events that will be produced
0075   std::vector<HepMC3::GenEvent> events;
0076   events.reserve(initialParticles.size());
0077
0078   for (const auto& part : initialParticles) {
0079     // Prepare the particle gun
0080     Geant4::HepMC3::PrimaryGeneratorAction::instance()->prepareParticleGun(
0081         part);
0082
0083     // Begin with the simulation
0084     m_runManager->BeamOn(1);
0085
0086     // Test if the event was aborted
0087     if (Geant4::HepMC3::SteppingAction::instance()->eventAborted()) {
0088       continue;
0089     }
0090
0091     // Set event start time
0092     HepMC3::GenEvent event = Geant4::HepMC3::EventAction::instance()->event();
0093     HepMC3::FourVector shift(0., 0., 0., part.time() / Acts::UnitConstants::mm);
0094     event.shift_position_by(shift);
0095
0096     // Set beam particle properties
0097     const Acts::Vector4 momentum4 =
0098         part.fourMomentum() / Acts::UnitConstants::GeV;
0099     HepMC3::FourVector beamMom4(momentum4[0], momentum4[1], momentum4[2],
0100                                 momentum4[3]);
0101     auto beamParticle = event.particles()[0];
0102     beamParticle->set_momentum(beamMom4);
0103     beamParticle->set_pid(part.pdg());
0104
0105     if (m_cfg.processSelect.empty()) {
0106       // Store the result
0107       events.push_back(std::move(event));
0108     } else {
0109       bool storeEvent = false;
0110       // Test if the event has a process of interest in it
0111       for (const auto& vertex : event.vertices()) {
0112         if (vertex->id() == -1) {
0113           vertex->add_particle_in(beamParticle);
0114         }
0115         const std::vector<std::string> vertexAttributes =
0116             vertex->attribute_names();
0117         for (const auto& att : vertexAttributes) {
0118           if ((vertex->attribute_as_string(att).find(m_cfg.processSelect) !=
0119                std::string::npos) &&
0120               !vertex->particles_in().empty() &&
0121               vertex->particles_in()[0]->attribute<HepMC3::IntAttribute>(
0122                   "TrackID") &&
0123               vertex->particles_in()[0]
0124                       ->attribute<HepMC3::IntAttribute>("TrackID")
0125                       ->value() == 1) {
0126             storeEvent = true;
0127             break;
0128           }
0129         }
0130         if (storeEvent) {
0131           break;
0132         }
0133       }
0134       // Store the result
0135       if (storeEvent) {
0136         // Remove vertices w/o outgoing particles and particles w/o production
0137         // vertices
0138         while (true) {
0139           bool sane = true;
0140           for (const auto& v : event.vertices()) {
0141             if (!v) {
0142               continue;
0143             }
0144             if (v->particles_out().empty()) {
0145               event.remove_vertex(v);
0146               sane = false;
0147             }
0148           }
0149           for (const auto& p : event.particles()) {
0150             if (!p) {
0151               continue;
0152             }
0153             if (!p->production_vertex()) {
0154               event.remove_particle(p);
0155               sane = false;
0156             }
0157           }
0158           if (sane) {
0159             break;
0160           }
0161         }
0162         events.push_back(std::move(event));
0163       }
0164     }
0165   }
0166
0167   ACTS_INFO(initialParticles.size() << " initial particles provided");
0168   ACTS_INFO(events.size() << " tracks generated");
0169
0170   // Write the recorded material to the event store
0171   m_outputEvents(context, std::move(events));
0172
0173   return ProcessCode::SUCCESS;
0174 }
0175
0176 }  // namespace ActsExamples