Geant4Decay.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Geant4Decay.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2021 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once

#include "Acts/Definitions/PdgParticle.hpp"
#include "Acts/Definitions/Units.hpp"
#include "ActsFatras/EventData/Particle.hpp"
#include "ActsFatras/Geant4/PDGtoG4Converter.hpp"

#include <cmath>
#include <limits>
#include <random>
#include <vector>

#include "G4ParticleDefinition.hh"

class G4RunManager;

namespace ActsFatras {

class Geant4Decay {
 public:
  using Scalar = Particle::Scalar;

  Geant4Decay();

  template <typename generator_t>
  Scalar generateProperTimeLimit(generator_t& generator,
                                 const Particle& particle) const;

  template <typename generator_t>
  std::vector<Particle> run(generator_t& generator, Particle& particle) const;

 private:
  std::vector<Particle> decayParticle(const Particle& parent) const;

  mutable G4RunManager* m_g4RunManager;  

  PDGtoG4Converter
      m_pdgToG4Conv;  
};

template <typename generator_t>
Particle::Scalar Geant4Decay::generateProperTimeLimit(
    generator_t& generator, const Particle& particle) const {
  // Get the particle properties
  const Acts::PdgParticle pdgCode = particle.pdg();
  // Keep muons stable
  if (makeAbsolutePdgParticle(pdgCode) == Acts::PdgParticle::eMuon) {
    return std::numeric_limits<Scalar>::infinity();
  }

  // Get the Geant4 particle
  G4ParticleDefinition* pDef = m_pdgToG4Conv.getParticleDefinition(pdgCode);

  // Fast exit if the particle is stable
  if (!pDef || pDef->GetPDGStable()) {
    return std::numeric_limits<Scalar>::infinity();
  }

  // Get average lifetime
  constexpr Scalar convertTime = Acts::UnitConstants::mm / CLHEP::s;
  const Scalar tau = pDef->GetPDGLifeTime() * convertTime;
  // Sample & return the lifetime
  std::uniform_real_distribution<Scalar> uniformDistribution{0., 1.};

  return -tau * std::log(uniformDistribution(generator));
}

template <typename generator_t>
std::vector<Particle> Geant4Decay::run(generator_t& /*generator*/,
                                       Particle& particle) const {
  return decayParticle(particle);
}
}  // namespace ActsFatras