Merge pull request key4hep#3 from delitez/propagatorAlg

EventGeneratorAlg

eventgen.py

@@ -0,0 +1,26 @@
+import os
+from pprint import pprint
+from Gaudi.Configuration import *
+
+from Configurables import EventGeneratorAlg
+
+sys.path.append('/home/delitez/ACTS/spack/k4actstracking')
+import actsUnits
+
+algList = []
+
+
+a = EventGeneratorAlg("MyEventGeneratorAlg")
+a.d0Sigma = 15 * actsUnits.um
+a.z0Sigma = 55 * actsUnits.mm
+a.tSigma = 1 * actsUnits.ns
+a.nMultiplicity = 5;
+a.nParticles = 10;
+a.objectPath = "/Event/MyParticle1"
+algList.append(a)
+
+
+
+from Configurables import ApplicationMgr
+
+ApplicationMgr(TopAlg=algList, EvtSel="NONE", EvtMax=1, ExtSvc=[], OutputLevel=DEBUG)

k4ActsTracking/src/components/EventGeneratorAlg.cpp

@@ -0,0 +1,139 @@
+#include "EventGeneratorAlg.h"
+#include "GaudiKernel/Service.h"
+
+
+using namespace Gaudi;
+
+DECLARE_COMPONENT(EventGeneratorAlg)
+
+EventGeneratorAlg::EventGeneratorAlg(const std::string& aName, ISvcLocator* aSvcLoc) : GaudiAlgorithm(aName, aSvcLoc) {}
+
+EventGeneratorAlg::~EventGeneratorAlg() { }
+
+StatusCode EventGeneratorAlg::initialize() {
+
+  return StatusCode::SUCCESS;
+}
+StatusCode EventGeneratorAlg::execute() {
+
+   SimParticleContainer particles;
+
+   static int seed = 5;
+   std::mt19937 rng{seed++};
+   std::normal_distribution<double> gauss(0., 1.);
+   std::round(gauss(rng));
+
+  size_t nPrimaryVertices = 0;
+
+  for (size_t n = nMultiplicity; 0 < n; --n) {
+    nPrimaryVertices+=1;
+
+// TODO : Random 4vec does not work. It is set to 0,0,0,0 temporarily.
+
+//    auto vertexPosition = (*vertex)(rng);
+Acts::Vector4 vertexPosition(0., 0., 0., 0.);
+
+
+    auto updateParticleInPlace = [&](ActsFatras::Particle& particle) {
+
+        const auto pid = ActsFatras::Barcode(particle.particleId())
+                                 .setVertexPrimary(nPrimaryVertices);
+
+        const auto pos4 = (vertexPosition + particle.fourPosition()).eval();
+
+        particle = particle.withParticleId(pid).setPosition4(pos4);
+    };
+
+
+     auto vertexParticles = genVertexParticles(rng,gauss);
+
+    for (auto& vertexParticle : vertexParticles) {
+      updateParticleInPlace(vertexParticle);
+    }
+
+    particles.merge(std::move(vertexParticles));
+
+  }
+
+
+  Container *containedParticle = new Container();
+  containedParticle->m_particles = particles;
+
+
+// This registiration must be revised because path is given in the python code and it crashes when it is not changed each time ---- TODO: Find a better way
+// Question: Should particles registered individually or all together  --- if all together, then the EventGeneratorAlg must run only once and it is okay as it is now.
+// Remark: Same name is okay if no compile, if compile, same name causes a crash
+
+  StatusCode sc = eventSvc()->registerObject(objectPath, containedParticle);
+
+  if( sc.isFailure() ) {
+     std::cout << "Object cannot be registered." << std::endl;
+     return StatusCode::FAILURE;
+  }
+  else{
+    std::cout << "Object is registered in " << objectPath << std::endl;
+  }
+
+  return StatusCode::SUCCESS;
+}
+StatusCode EventGeneratorAlg::finalize() {
+
+  return StatusCode::SUCCESS;
+}
+
+
+SimParticleContainer EventGeneratorAlg::genVertexParticles(std::mt19937& rng, std::normal_distribution<double>& gauss) {
+
+  using UniformIndex = std::uniform_int_distribution<unsigned int>;
+
+
+  std::uniform_real_distribution<double> phiDist(0, 2 * M_PI);
+  std::uniform_real_distribution<double> etaDist(-4.0, 4.0);
+  std::uniform_real_distribution<double> ptDist(10., 20.);
+  std::uniform_real_distribution<double> qDist(0., 1.);
+
+  double d0     = d0Sigma * gauss(rng);
+  double z0     = z0Sigma * gauss(rng);
+  double charge = qDist(rng) > 0.5 ? 1. : -1.;
+  double t      = tSigma * gauss(rng);
+
+  UniformIndex particleTypeChoice(0u, qDist(rng) ? 1u : 0u);
+  Acts::PdgParticle pdg = Acts::PdgParticle::eMuon;
+
+  const Acts::PdgParticle pdgChoices[] = {pdg, static_cast<Acts::PdgParticle>(-pdg), };    //warum leer zeichen?
+                                                                                          // braucht man 3 elements statt 2?
+  const double qChoices[] = {charge,-charge, };
+
+  SimParticleContainer particles;
+
+  for (size_t ip = 1 ; ip <= nParticles; ++ip) {
+    const auto pid = ActsFatras::Barcode(0u).setParticle(ip);
+    const unsigned int type = particleTypeChoice(rng);
+    const double q = qChoices[type];
+    const double phi    = phiDist(rng);
+    double eta    = etaDist(rng);
+    double theta  = 2 * atan(exp(-eta));
+    double pt     = ptDist(rng);
+    double p      = pt / sin(theta);
+
+    Acts::Vector3 direction;
+
+    direction[Acts::eMom0] = sin(theta) * cos(phi);
+    direction[Acts::eMom1] = sin(theta) * sin(phi);
+    direction[Acts::eMom2] = cos(theta);
+
+    //TODO: Fix the mass -- how to extract it from pdg? config file? by hand?
+
+    double mass = 0.5;
+
+    ActsFatras::Particle particle(pid, pdg, q, mass);
+    particle.setDirection(direction);
+    particle.setAbsoluteMomentum(p);
+
+    particles.insert(particles.end(), std::move(particle));
+
+  }
+
+return particles;
+
+};

k4ActsTracking/src/components/EventGeneratorAlg.h

@@ -0,0 +1,130 @@
+#ifndef EventGeneratorSvc_H
+#define EventGeneratorSvc_H
+
+
+#include "Acts/Utilities/Helpers.hpp"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/Service.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "Gaudi/Property.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "IEventGeneratorSvc.h"
+#include "Acts/Definitions/Algebra.hpp"
+#include "ActsFatras/EventData/Particle.hpp"
+#include "ActsFatras/Utilities/ParticleData.hpp"
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/Utilities/PdgParticle.hpp"
+#include <boost/container/flat_set.hpp>
+#include <random>
+#include <cmath>
+
+namespace detail {
+struct CompareParticleId {
+  using is_transparent = void;
+  constexpr bool operator()(const ActsFatras::Particle& lhs,
+                            const ActsFatras::Particle& rhs) const {
+    return lhs.particleId() < rhs.particleId();
+  }
+  constexpr bool operator()(ActsFatras::Barcode lhs,
+                            const ActsFatras::Particle& rhs) const {
+    return lhs < rhs.particleId();
+  }
+  constexpr bool operator()(const ActsFatras::Particle& lhs,
+                            ActsFatras::Barcode rhs) const {
+    return lhs.particleId() < rhs;
+  }
+};
+}
+
+using SimParticleContainer =
+    ::boost::container::flat_set<::ActsFatras::Particle,
+                                 detail::CompareParticleId>;
+
+
+
+class EventGeneratorAlg : public GaudiAlgorithm {
+
+public:
+
+SimParticleContainer genVertexParticles(std::mt19937& rng, std::normal_distribution<double>& gauss);
+
+
+  struct MultiplicityGenerator {
+    virtual ~MultiplicityGenerator() = default;
+    virtual size_t operator()(std::mt19937& rng) const = 0;
+  };
+
+  struct VertexGenerator{
+    virtual ~VertexGenerator() = default;
+    virtual Acts::Vector4 operator()(std::mt19937& rng) const = 0;
+  };
+
+  struct Container  : public DataObject {
+    SimParticleContainer m_particles;
+  };
+
+
+private:
+
+
+public:
+  Gaudi::Property<double> d0Sigma{this, "d0Sigma", 0, "Option for d0 gaussian sigma"};
+
+  Gaudi::Property<double> z0Sigma{this, "z0Sigma", 0, "Option for z0 gaussian sigma"};
+
+  Gaudi::Property<double> phiSigma{this, "phiSigma", 0,
+                                   "Option for phi gaussian sigma (used for covariance transport)"};
+  Gaudi::Property<double> tSigma{this, "tSigma", 0, "Option for t gaussian sigma (used for covariance transport)"};
+
+  Gaudi::Property<int> nParticles{this, "nParticles", 0, "Number of particles."};
+
+  Gaudi::Property<int> nMultiplicity{this, "nMultiplicity", 0, "Number of multiplicity."};
+
+  Gaudi::Property<std::string> objectPath{this, "objectPath", " ", "Path for the object."};
+
+  std::shared_ptr<VertexGenerator> vertex = nullptr;
+
+  EventGeneratorAlg(const std::string& name, ISvcLocator* svc);
+
+  virtual ~EventGeneratorAlg();
+
+  virtual StatusCode initialize() final;
+
+  virtual StatusCode execute() final;
+
+  virtual StatusCode finalize() final;
+
+
+};
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#endif

k4ActsTracking/src/components/objectTEST.cpp

@@ -0,0 +1,44 @@
+#include "objectTEST.h"
+
+#include <iostream>
+#include <string>
+#include <filesystem>
+#include <unistd.h>
+#include "PropagatorAlg.h"
+
+//using std::filesystem::current_path;
+
+
+DECLARE_COMPONENT(objectTest)
+
+objectTest::objectTest(const std::string& aName, ISvcLocator* aSvcLoc) : GaudiAlgorithm(aName, aSvcLoc) {}
+
+objectTest::~objectTest() {}
+
+StatusCode objectTest::initialize() {return StatusCode::SUCCESS; }
+
+StatusCode objectTest::execute() {
+
+  typedef std::vector<int> MyTestVector;
+  DataObject *pObject;
+  //std::string objectPath = "./testDir";
+
+  StatusCode sc = eventSvc()->registerObject(".", pObject);
+
+   if( sc.isFailure() ) {
+      std::cout << "CANNOT initialize sc" << std::endl;
+      return StatusCode::FAILURE;
+   }
+   else{
+     std::cout << "CAN initialize sc" << std::endl;
+   }
+
+   MyTestVector *tv = 0;
+   tv = dynamic_cast<MyTestVector *> (pObject);
+
+
+   std::cout << "Object Test is alive!" << std::endl;
+  return StatusCode::SUCCESS;
+}
+
+StatusCode objectTest::finalize() { return StatusCode::SUCCESS; }

k4ActsTracking/src/components/objectTEST.h

@@ -0,0 +1,28 @@
+#pragma once
+
+// GAUDI
+#include "Gaudi/Property.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/IDataProviderSvc.h"
+//#include "GaudiKernel/RegistryEntry.h"
+
+class objectTest : public GaudiAlgorithm {
+public:
+  explicit objectTest(const std::string&, ISvcLocator*);
+  virtual ~objectTest();
+  /**  Initialize.
+   *   @return status code
+   */
+  virtual StatusCode initialize() final;
+  /**  Execute.
+   *   @return status code
+   */
+  virtual StatusCode execute() final;
+  /**  Finalize.
+   *   @return status code
+   */
+  virtual StatusCode finalize() final;
+
+private:
+
+};