MC2DummyParticle.cpp 4.46 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

#include "JugBase/DataHandle.h"
10
#include "JugBase/UniqueID.h"
11
12
13
14

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

17
namespace Jug::Base {
18

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

29
30
      const int32_t kMonteCarloSource{uniqueID<int32_t>("mcparticles")};

31
32
      MC2DummyParticle(const std::string& name, ISvcLocator* svcLoc) 
        : GaudiAlgorithm(name, svcLoc)
33
        , AlgorithmIDMixin(name, info())
34
35
36
      {
        declareProperty("inputCollection", m_inputHitCollection, "mcparticles");
        declareProperty("outputCollection", m_outputHitCollection, "DummyReconstructedParticles");
37
        declareProperty("outputRelations", m_outputRelCollection, "DummyReconstructedParticles");
38
39
40
41
42
      }
      StatusCode initialize() override
      {
        if (GaudiAlgorithm::initialize().isFailure())
          return StatusCode::FAILURE;
Whitney Armstrong's avatar
Whitney Armstrong committed
43
44

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

          // 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.
Sylvester Joosten's avatar
Sylvester Joosten committed
68
69
70
71
72
73
74
75
76
77
          const double pgen     = p.ps().mag();
          const float momentum = pgen * m_gaussDist();
          const float energy   = p.energy();
          const float px       = p.ps().x * momentum / pgen;
          const float py       = p.ps().y * momentum / pgen;
          const float pz       = p.ps().z * momentum / pgen;
          // @TODO: vertex smearing
          const float vx       = p.vs().x;
          const float vy       = p.vs().y;
          const float vz       = p.vs().z;
78

79
80
          eic::VectorXYZ psmear{px, py, pz};

Sylvester Joosten's avatar
Sylvester Joosten committed
81
          eic::ReconstructedParticle rec_part{
82
            {ID++, algorithmID()},            // Unique index
83
            psmear,                           // 3-momentum [GeV]
Sylvester Joosten's avatar
Sylvester Joosten committed
84
            {vx, vy, vz},                     // Vertex [mm]
85
            static_cast<float>(p.time()),     // time [ns]
Sylvester Joosten's avatar
Sylvester Joosten committed
86
87
88
            p.pdgID(),                        // PDG type
            static_cast<int16_t>(p.status()), // Status
            static_cast<int16_t>(p.charge()), // Charge
89
            1.,                               // particle weight
90
            {psmear.theta(), psmear.phi()},   // direction
Sylvester Joosten's avatar
Sylvester Joosten committed
91
92
            momentum,                         // 3-momentum magnitude [GeV]
            energy,                           // energy [GeV]
93
            static_cast<float>(p.mass())};    // mass [GeV]
Sylvester Joosten's avatar
Sylvester Joosten committed
94

95
          out_parts->push_back(rec_part);
96
97
98
99

          eic::ReconstructedParticleRelations rel;
          rel.mcID({p.ID(), kMonteCarloSource});
          relations->push_back(rel);
100
101
102
103
104
105
106
        }
        return StatusCode::SUCCESS;
      }
    };

    DECLARE_COMPONENT(MC2DummyParticle)

107
} // namespace Jug::Base
108