MC2DummyParticle.cpp 3.43 KB
Newer Older
1
#include "Gaudi/Algorithm.h"
2
3
4
#include "GaudiAlg/GaudiTool.h"
#include "GaudiAlg/Producer.h"
#include "GaudiAlg/Transformer.h"
Whitney Armstrong's avatar
Whitney Armstrong committed
5
#include "GaudiKernel/RndmGenerators.h"
6
7
#include <algorithm>
#include <cmath>
8
9
10
11
12
13
14
15
16
17
18
19
20

// FCCSW
#include "JugBase/DataHandle.h"

// Event Model related classes
#include "dd4pod/Geant4ParticleCollection.h"
#include "eicd/ReconstructedParticleCollection.h"

namespace Jug {
  namespace Base {

    class MC2DummyParticle : public GaudiAlgorithm {
    public:
21
22
23
24
25
      DataHandle<dd4pod::Geant4ParticleCollection> m_inputHitCollection{"mcparticles", Gaudi::DataHandle::Reader, this};
      DataHandle<eic::ReconstructedParticleCollection> m_outputHitCollection{"DummyReconstructedParticles",
                                                                             Gaudi::DataHandle::Writer, this};
      Rndm::Numbers                                    m_gaussDist;
      Gaudi::Property<double>                          m_smearing{this, "smearing", 0.01 /* 1 percent*/};
26

27
28
29
30
31
32
33
34
35
      MC2DummyParticle(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
      {
        declareProperty("inputCollection", m_inputHitCollection, "mcparticles");
        declareProperty("outputCollection", m_outputHitCollection, "DummyReconstructedParticles");
      }
      StatusCode initialize() override
      {
        if (GaudiAlgorithm::initialize().isFailure())
          return StatusCode::FAILURE;
Whitney Armstrong's avatar
Whitney Armstrong committed
36
37

        IRndmGenSvc* randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
Whitney Armstrong's avatar
Whitney Armstrong committed
38
        StatusCode   sc      = m_gaussDist.initialize(randSvc, Rndm::Gauss(1.0, m_smearing.value()));
Whitney Armstrong's avatar
Whitney Armstrong committed
39
40
41
        if (!sc.isSuccess()) {
          return StatusCode::FAILURE;
        }
42
43
44
45
46
47
48
49
        return StatusCode::SUCCESS;
      }
      StatusCode execute() override
      {
        // input collection
        const dd4pod::Geant4ParticleCollection* parts = m_inputHitCollection.get();
        // output collection
        auto out_parts = m_outputHitCollection.createAndPut();
Sylvester Joosten's avatar
Sylvester Joosten committed
50
        int ID = 0;
51
        for (const auto& p : *parts) {
52
53
54
          if (p.genStatus() != 1) {
            continue;
          }
55
56
57
58
59
60
61
62
63
64
65

          // for now just use total momentum smearing as this is the largest effect,
          // ideally we should also smear the angles but this should be good enough
          // for now.
          const double pgen     = std::hypot(p.psx(), p.psy(), p.psz());
          const double momentum = pgen * m_gaussDist();
          const double energy   = std::hypot(momentum, p.mass());
          const double px       = p.psx() * momentum / pgen;
          const double py       = p.psy() * momentum / pgen;
          const double pz       = p.psz() * momentum / pgen;

Sylvester Joosten's avatar
Sylvester Joosten committed
66
67
68
69
70
71
72
73
74
75
76
77
78
          eic::ReconstructedParticle rec_part{
            ID++,                     // Unique index
            {px, py, pz},             // 3-momentum [GeV]
            {0, 0, 0},                // @TODO: Vertex [mm]
            0.,                       // @TODO: time [ns]
            p.pdgID(),                // PDG type
            static_cast<int16_t>(0),  // @TODO: Status
            p.charge(),               // Charge
            momentum,                 // 3-momentum magnitude [GeV]
            energy,                   // energy [GeV]
            p.mass(),                 // mass [GeV]
            1.};                      // particle weight

79
          out_parts->push_back(rec_part);
80
81
82
83
84
85
86
87
        }
        return StatusCode::SUCCESS;
      }
    };

    DECLARE_COMPONENT(MC2DummyParticle)

  } // namespace Base
88
} // namespace Jug
89