CalorimeterHitDigi.cpp 9.55 KB
Newer Older
Chao Peng's avatar
Chao Peng committed
1
2
3
4
// A general digitization for CalorimeterHit from simulation
// 1. Smear energy deposit with a/sqrt(E/GeV) + b + c/E or a/sqrt(E/GeV) (relative value)
// 2. Digitize the energy with dynamic ADC range and add pedestal (mean +- sigma)
// 3. Time conversion with smearing resolution (absolute value)
5
// 4. Signal is summed if the SumFields are provided
Chao Peng's avatar
Chao Peng committed
6
7
8
9
10
11
//
// Author: Chao Peng
// Date: 06/02/2021

#include <algorithm>
#include <cmath>
12
#include <unordered_map>
Chao Peng's avatar
Chao Peng committed
13
14
15
16
17
18
19

#include "GaudiAlg/GaudiTool.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "Gaudi/Property.h"
#include "GaudiKernel/RndmGenerators.h"

20
21
22
23
#include "DDRec/CellIDPositionConverter.h"
#include "DDSegmentation/BitFieldCoder.h"

#include "JugBase/IGeoSvc.h"
Chao Peng's avatar
Chao Peng committed
24
#include "JugBase/DataHandle.h"
25
#include "JugBase/UniqueID.h"
Chao Peng's avatar
Chao Peng committed
26

27
28
29
#include "fmt/format.h"
#include "fmt/ranges.h"

Chao Peng's avatar
Chao Peng committed
30
31
32
33
34
// Event Model related classes
#include "dd4pod/CalorimeterHitCollection.h"
#include "eicd/RawCalorimeterHitCollection.h"
#include "eicd/RawCalorimeterHitData.h"

35

Chao Peng's avatar
Chao Peng committed
36
37
38
39
using namespace Gaudi::Units;

namespace Jug::Digi {

Whitney Armstrong's avatar
Whitney Armstrong committed
40
41
42
43
44
  /** Generic calorimeter hit digitiziation.
   *
   * \ingroup digi
   * \ingroup calorimetry
   */
45
  class CalorimeterHitDigi : public GaudiAlgorithm, AlgorithmIDMixin<> {
Whitney Armstrong's avatar
Whitney Armstrong committed
46
  public:
Chao Peng's avatar
Chao Peng committed
47
    // additional smearing resolutions
Whitney Armstrong's avatar
Whitney Armstrong committed
48
    Gaudi::Property<std::vector<double>> u_eRes{this, "energyResolutions", {}}; // a/sqrt(E/GeV) + b + c/(E/GeV)
49
    Gaudi::Property<double>              m_tRes{this, "timeResolution", 0.0 * ns};
50

Chao Peng's avatar
Chao Peng committed
51
    // digitization settings
52
53
54
55
56
57
    Gaudi::Property<int>                m_capADC{this, "capacityADC", 8096};
    Gaudi::Property<double>             m_dyRangeADC{this, "dynamicRangeADC", 100 * MeV};
    Gaudi::Property<int>                m_pedMeanADC{this, "pedestalMean", 400};
    Gaudi::Property<double>             m_pedSigmaADC{this, "pedestalSigma", 3.2};
    Gaudi::Property<double>             m_resolutionTDC{this, "resolutionTDC", 10 * ps};

58
59
60
61
62
63
    Gaudi::Property<double>             m_corrMeanScale{this, "scaleResponse", 1.0};
    // These are better variable names for the "energyResolutions" array which is a bit
    // magic @FIXME
    //Gaudi::Property<double>             m_corrSigmaCoeffE{this, "responseCorrectionSigmaCoeffE", 0.0};
    //Gaudi::Property<double>             m_corrSigmaCoeffSqrtE{this, "responseCorrectionSigmaCoeffSqrtE", 0.0};

64
65
66
67
68
69
70
71
72
    // signal sums
    // @TODO: implement signal sums with timing
    // field names to generate id mask, the hits will be grouped by masking the field
    Gaudi::Property<std::vector<std::string>> u_fields{this, "signalSumFields", {}};
    // ref field ids are used for the merged hits, 0 is used if nothing provided
    Gaudi::Property<std::vector<int>>         u_refs{this, "fieldRefNumbers", {}};
    Gaudi::Property<std::string>              m_geoSvcName{this, "geoServiceName", "GeoSvc"};
    Gaudi::Property<std::string>              m_readout{this, "readoutClass", ""};

73
    // unitless counterparts of inputs
74
75
76
77
78
79
80
81
82
    double           dyRangeADC, stepTDC, tRes, eRes[3] = {0., 0., 0.};
    Rndm::Numbers    m_normDist;
    SmartIF<IGeoSvc> m_geoSvc;
    uint64_t         id_mask, ref_mask;

    DataHandle<dd4pod::CalorimeterHitCollection> m_inputHitCollection{
      "inputHitCollection", Gaudi::DataHandle::Reader, this};
    DataHandle<eic::RawCalorimeterHitCollection> m_outputHitCollection{
      "outputHitCollection", Gaudi::DataHandle::Writer, this};
Chao Peng's avatar
Chao Peng committed
83

Whitney Armstrong's avatar
Whitney Armstrong committed
84
    //  ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
85
    CalorimeterHitDigi(const std::string& name, ISvcLocator* svcLoc)
86
87
      : GaudiAlgorithm(name, svcLoc)
      , AlgorithmIDMixin{name, info()}
Chao Peng's avatar
Chao Peng committed
88
    {
Whitney Armstrong's avatar
Whitney Armstrong committed
89
90
      declareProperty("inputHitCollection", m_inputHitCollection, "");
      declareProperty("outputHitCollection", m_outputHitCollection, "");
Chao Peng's avatar
Chao Peng committed
91
92
93
94
    }

    StatusCode initialize() override
    {
Whitney Armstrong's avatar
Whitney Armstrong committed
95
96
97
98
99
100
101
102
103
104
105
106
107
      if (GaudiAlgorithm::initialize().isFailure()) {
        return StatusCode::FAILURE;
      }
      // random number generator from service
      auto randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
      auto sc      = m_normDist.initialize(randSvc, Rndm::Gauss(0.0, 1.0));
      if (!sc.isSuccess()) {
        return StatusCode::FAILURE;
      }
      // set energy resolution numbers
      for (size_t i = 0; i < u_eRes.size() && i < 3; ++i) {
        eRes[i] = u_eRes[i];
      }
108

Whitney Armstrong's avatar
Whitney Armstrong committed
109
110
111
      // using juggler internal units (GeV, mm, radian, ns)
      dyRangeADC = m_dyRangeADC.value() / GeV;
      tRes       = m_tRes.value() / ns;
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
      stepTDC    = ns / m_resolutionTDC.value();

      // need signal sum
      if (u_fields.value().size()) {
        m_geoSvc = service(m_geoSvcName);
        // sanity checks
        if (!m_geoSvc) {
          error() << "Unable to locate Geometry Service. "
                  << "Make sure you have GeoSvc and SimSvc in the right order in the configuration."
                  << endmsg;
          return StatusCode::FAILURE;
        }
        if (m_readout.value().empty()) {
          error() << "readoutClass is not provided, it is needed to know the fields in readout ids"
                  << endmsg;
          return StatusCode::FAILURE;
        }

        // get decoders
        try {
          auto id_desc = m_geoSvc->detector()->readout(m_readout).idSpec();
          id_mask      = 0;
          std::vector<std::pair<std::string, int>> ref_fields;
          for (size_t i = 0; i < u_fields.size(); ++i) {
            id_mask |= id_desc.field(u_fields[i])->mask();
            // use the provided id number to find ref cell, or use 0
            int ref = i < u_refs.size() ? u_refs[i] : 0;
            ref_fields.push_back({u_fields[i], ref});
          }
          ref_mask = id_desc.encode(ref_fields);
          // debug() << fmt::format("Referece id mask for the fields {:#064b}", ref_mask) << endmsg;
        } catch (...) {
          error() << "Failed to load ID decoder for " << m_readout << endmsg;
          return StatusCode::FAILURE;
        }
        id_mask = ~id_mask;
        info() << fmt::format("ID mask in {:s}: {:#064b}", m_readout, id_mask) << endmsg;
        return StatusCode::SUCCESS;
      }
151

Whitney Armstrong's avatar
Whitney Armstrong committed
152
      return StatusCode::SUCCESS;
Chao Peng's avatar
Chao Peng committed
153
154
155
156
    }

    StatusCode execute() override
    {
157
158
159
160
161
162
163
164
165
166
      if (u_fields.value().size()) {
        signal_sum_digi();
      } else {
        single_hits_digi();
      }
      return StatusCode::SUCCESS;
    }

  private:
    void single_hits_digi() {
Whitney Armstrong's avatar
Whitney Armstrong committed
167
168
169
170
      // input collections
      const auto simhits = m_inputHitCollection.get();
      // Create output collections
      auto rawhits = m_outputHitCollection.createAndPut();
Sylvester Joosten's avatar
Sylvester Joosten committed
171
      int nhits = 0;
Whitney Armstrong's avatar
Whitney Armstrong committed
172
173
      for (const auto& ahit : *simhits) {
        // Note: juggler internal unit of energy is GeV
174
175
176
177
178
179
180
181
182
183
184
185
186
187
        // --> optional energy correction to allow for modifying the response for
        // effective impelentations, e.g. a homogenous implementation of a fiber
        // calorimeter.
        const double eDep    = ahit.energyDeposit() * m_corrMeanScale;

        // apply additional calorimeter noise to corrected energy deposit
        const double eResRel = (eDep > 1e-6)
                                   ? m_normDist() * std::sqrt(std::pow(eRes[0] / std::sqrt(eDep), 2) +
                                                              std::pow(eRes[1], 2) + std::pow(eRes[2] / (eDep), 2))
                                   : 0;

        const double ped    = m_pedMeanADC + m_normDist() * m_pedSigmaADC;
        const long long adc = std::llround(ped + eDep * (1. + eResRel) / dyRangeADC * m_capADC);
        const long long tdc = std::llround((ahit.truth().time + m_normDist() * tRes) * stepTDC);
188

Sylvester Joosten's avatar
Sylvester Joosten committed
189
        eic::RawCalorimeterHit rawhit(
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
          {nhits++, algorithmID()},
          (long long)ahit.cellID(),
          (adc > m_capADC.value() ? m_capADC.value() : adc),
          tdc
        );
        rawhits->push_back(rawhit);
      }
    }

    void signal_sum_digi() {
      const auto simhits = m_inputHitCollection.get();
      auto rawhits = m_outputHitCollection.createAndPut();

      // find the hits that belong to the same group (for merging)
      std::unordered_map<long long, std::vector<dd4pod::ConstCalorimeterHit>> merge_map;
      for (const auto &ahit : *simhits) {
        int64_t hid = (ahit.cellID() & id_mask) | ref_mask;
        auto    it  = merge_map.find(hid);

        if (it == merge_map.end()) {
          merge_map[hid] = {ahit};
        } else {
          it->second.push_back(ahit);
        }
      }

      // signal sum
      int nhits = 0;
      for (auto &[id, hits] : merge_map) {
        double edep     = hits[0].energyDeposit();
        double time     = hits[0].truth().time;
        double max_edep = hits[0].energyDeposit();
        // sum energy, take time from the most energetic hit
        for (size_t i = 1; i < hits.size(); ++i) {
          edep += hits[i].energyDeposit();
          if (hits[i].energyDeposit() > max_edep) {
            max_edep = hits[i].energyDeposit();
            time     = hits[i].truth().time;
          }
        }
230
231
232
233
234
235
236
237

        double eResRel = 0.;
        // safety check
        if (edep > 1e-6) {
            eResRel = m_normDist() * eRes[0] / std::sqrt(edep) +
                      m_normDist() * eRes[1] +
                      m_normDist() * eRes[2] / edep;
        }
238
239
240
241
242
243
244
245
246
247
        double    ped     = m_pedMeanADC + m_normDist() * m_pedSigmaADC;
        long long adc     = std::llround(ped + edep * (1. + eResRel) / dyRangeADC * m_capADC);
        long long tdc     = std::llround((time + m_normDist() * tRes) * stepTDC);

        eic::RawCalorimeterHit rawhit(
          {nhits++, algorithmID()},
          id,
          (adc > m_capADC.value() ? m_capADC.value() : adc),
          tdc
        );
Whitney Armstrong's avatar
Whitney Armstrong committed
248
249
        rawhits->push_back(rawhit);
      }
Chao Peng's avatar
Chao Peng committed
250
    }
Whitney Armstrong's avatar
Whitney Armstrong committed
251
252
  };
  DECLARE_COMPONENT(CalorimeterHitDigi)
Chao Peng's avatar
Chao Peng committed
253
254

} // namespace Jug::Digi