CalorimeterHitReco.cpp 8.91 KB
Newer Older
1
2
3
4
// Reconstruct digitized outputs, paired with Jug::Digi::CalorimeterHitDigi
// Author: Chao Peng
// Date: 06/14/2021

5
6
#include "fmt/format.h"
#include "fmt/ranges.h"
Sylvester Joosten's avatar
Sylvester Joosten committed
7
8
#include <algorithm>
#include <bitset>
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

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

#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/Surface.h"
#include "DDRec/SurfaceManager.h"

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

// Event Model related classes
#include "eicd/CalorimeterHitCollection.h"
#include "eicd/RawCalorimeterHitCollection.h"

using namespace Gaudi::Units;

namespace Jug::Reco {

Whitney Armstrong's avatar
Whitney Armstrong committed
34
35
36
37
38
39

  /** Calorimeter hit reconstruction.
   *
   * Reconstruct digitized outputs, paired with Jug::Digi::CalorimeterHitDigi
   * \ingroup reco
   */
Sylvester Joosten's avatar
Sylvester Joosten committed
40
41
  class CalorimeterHitReco : public GaudiAlgorithm {
  public:
42
43
44
45
46
    // length unit from dd4hep, should be fixed
    Gaudi::Property<double> m_lUnit{this, "lengthUnit", dd4hep::mm};

    // digitization settings, must be consistent with digi class
    Gaudi::Property<int>    m_capADC{this, "capacityADC", 8096};
Sylvester Joosten's avatar
Sylvester Joosten committed
47
    Gaudi::Property<double> m_dyRangeADC{this, "dynamicRangeADC", 100. * MeV};
48
49
50
51
52
    Gaudi::Property<int>    m_pedMeanADC{this, "pedestalMean", 400};
    Gaudi::Property<double> m_pedSigmaADC{this, "pedestalSigma", 3.2};

    // zero suppression values
    Gaudi::Property<double> m_thresholdADC{this, "thresholdFactor", 3.0};
53

Sylvester Joosten's avatar
Sylvester Joosten committed
54
55
56
    // Calibration!
    Gaudi::Property<double> m_sampFrac{this, "samplingFraction", 1.0};

57
58
    // unitless counterparts of the input parameters
    double                  dyRangeADC;
59

Sylvester Joosten's avatar
Sylvester Joosten committed
60
61
62
63
    DataHandle<eic::RawCalorimeterHitCollection> m_inputHitCollection{
        "inputHitCollection", Gaudi::DataHandle::Reader, this};
    DataHandle<eic::CalorimeterHitCollection> m_outputHitCollection{
        "outputHitCollection", Gaudi::DataHandle::Writer, this};
64
65

    // geometry service to get ids, ignored if no names provided
Sylvester Joosten's avatar
Sylvester Joosten committed
66
67
68
69
70
71
72
    Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
    Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};
    Gaudi::Property<std::string> m_layerField{this, "layerField", ""};
    Gaudi::Property<std::string> m_sectorField{this, "sectorField", ""};
    SmartIF<IGeoSvc>             m_geoSvc;
    dd4hep::BitFieldCoder*       id_dec = nullptr;
    size_t                       sector_idx, layer_idx;
73

74
75
    // name of detelment or fields to find the local detector (for global->local transform)
    // if nothing is provided, the lowest level DetElement (from cellID) will be used
Sylvester Joosten's avatar
Sylvester Joosten committed
76
    Gaudi::Property<std::string>              m_localDetElement{this, "localDetElement", ""};
77
    Gaudi::Property<std::vector<std::string>> u_localDetFields{this, "localDetFields", {}};
Chao Peng's avatar
Chao Peng committed
78
    dd4hep::DetElement                        local;
Sylvester Joosten's avatar
Sylvester Joosten committed
79
    size_t                                    local_mask = ~0;
80

81
82
    CalorimeterHitReco(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
    {
Sylvester Joosten's avatar
Sylvester Joosten committed
83
84
      declareProperty("inputHitCollection", m_inputHitCollection, "");
      declareProperty("outputHitCollection", m_outputHitCollection, "");
85
86
87
88
    }

    StatusCode initialize() override
    {
Sylvester Joosten's avatar
Sylvester Joosten committed
89
90
91
92
93
94
95
96
97
98
99
100
      if (GaudiAlgorithm::initialize().isFailure()) {
        return StatusCode::FAILURE;
      }

      m_geoSvc = service(m_geoSvcName);
      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;
      }

101
102
103
      // unitless conversion
      dyRangeADC = m_dyRangeADC.value()/GeV;

Sylvester Joosten's avatar
Sylvester Joosten committed
104
105
106
107
      // do not get the layer/sector ID if no readout class provided
      if (m_readout.value().empty()) {
        return StatusCode::SUCCESS;
      }
108

Sylvester Joosten's avatar
Sylvester Joosten committed
109
110
111
112
113
      auto id_spec = m_geoSvc->detector()->readout(m_readout).idSpec();
      try {
        id_dec = id_spec.decoder();
        if (m_sectorField.value().size()) {
          sector_idx = id_dec->index(m_sectorField);
114
          info() << "Find sector field " << m_sectorField.value() << ", index = " << sector_idx << endmsg;
115
        }
Sylvester Joosten's avatar
Sylvester Joosten committed
116
117
        if (m_layerField.value().size()) {
          layer_idx = id_dec->index(m_layerField);
118
          info() << "Find layer field " << m_layerField.value() << ", index = " << sector_idx << endmsg;
119
        }
Sylvester Joosten's avatar
Sylvester Joosten committed
120
121
122
123
      } catch (...) {
        error() << "Failed to load ID decoder for " << m_readout << endmsg;
        return StatusCode::FAILURE;
      }
124

Sylvester Joosten's avatar
Sylvester Joosten committed
125
126
      // local detector name has higher priority
      if (m_localDetElement.value().size()) {
127
        try {
Chao Peng's avatar
Chao Peng committed
128
          local = m_geoSvc->detector()->detector(m_localDetElement.value());
Sylvester Joosten's avatar
Sylvester Joosten committed
129
130
          info() << "Local coordinate system from DetElement " << m_localDetElement.value()
                 << endmsg;
131
        } catch (...) {
Sylvester Joosten's avatar
Sylvester Joosten committed
132
133
134
          error() << "Failed to locate local coordinate system from DetElement "
                  << m_localDetElement.value() << endmsg;
          return StatusCode::FAILURE;
135
        }
136
      // or get from fields
Sylvester Joosten's avatar
Sylvester Joosten committed
137
138
139
140
      } else {
        std::vector<std::pair<std::string, int>> fields;
        for (auto& f : u_localDetFields.value()) {
          fields.push_back({f, 0});
141
        }
Sylvester Joosten's avatar
Sylvester Joosten committed
142
143
144
145
146
147
148
149
150
        local_mask = id_spec.get_mask(fields);
        // use all fields if nothing provided
        if (fields.empty()) {
          local_mask = ~0;
        }
        info() << fmt::format("Local DetElement mask {:#064b} from fields [{}]", local_mask,
                              fmt::join(fields, ", "))
               << endmsg;
      }
151

Sylvester Joosten's avatar
Sylvester Joosten committed
152
      return StatusCode::SUCCESS;
153
154
155
156
    }

    StatusCode execute() override
    {
Sylvester Joosten's avatar
Sylvester Joosten committed
157
158
159
160
      // input collections
      const auto& rawhits = *m_inputHitCollection.get();
      // create output collections
      auto& hits = *m_outputHitCollection.createAndPut();
Chao Peng's avatar
Chao Peng committed
161
      auto converter = m_geoSvc->cellIDPositionConverter();
Sylvester Joosten's avatar
Sylvester Joosten committed
162
163
164
165
166
167

      // energy time reconstruction
      for (const auto& rh : rawhits) {
        // did not pass the zero-suppression threshold
        if ((rh.amplitude() - m_pedMeanADC) < m_thresholdADC * m_pedSigmaADC) {
          continue;
168
169
        }

Sylvester Joosten's avatar
Sylvester Joosten committed
170
        // convert ADC -> energy
Sylvester Joosten's avatar
Sylvester Joosten committed
171
        float energy = (rh.amplitude() - m_pedMeanADC) / static_cast<float>(m_capADC.value()) * dyRangeADC / m_sampFrac;
Sylvester Joosten's avatar
Sylvester Joosten committed
172

Sylvester Joosten's avatar
Sylvester Joosten committed
173
174
        float time = rh.time()*1.e-6; // ns @FIXME: this should not be hardcoded
        auto  cellID   = rh.cellID();
175
        int   lid  = ((id_dec != nullptr) && m_layerField.value().size())
Sylvester Joosten's avatar
Sylvester Joosten committed
176
                      ? static_cast<int>(id_dec->get(cellID, layer_idx))
Sylvester Joosten's avatar
Sylvester Joosten committed
177
                      : -1;
178
        int   sid  = ((id_dec != nullptr) && m_sectorField.value().size())
Sylvester Joosten's avatar
Sylvester Joosten committed
179
                      ? static_cast<int>(id_dec->get(cellID, sector_idx))
Sylvester Joosten's avatar
Sylvester Joosten committed
180
181
                      : -1;
        // global positions
Sylvester Joosten's avatar
Sylvester Joosten committed
182
        auto gpos = converter->position(cellID);
Sylvester Joosten's avatar
Sylvester Joosten committed
183
184
185
        // local positions
        if (m_localDetElement.value().empty()) {
          auto volman = m_geoSvc->detector()->volumeManager();
Sylvester Joosten's avatar
Sylvester Joosten committed
186
          local       = volman.lookupDetElement(cellID & local_mask);
Sylvester Joosten's avatar
Sylvester Joosten committed
187
        }
Chao Peng's avatar
Chao Peng committed
188
        auto pos = local.nominal().worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
189
        // cell dimension
Chao Peng's avatar
Chao Peng committed
190
191
192
        std::vector<double> cdim;
        // get segmentation dimensions
        if (converter->findReadout(local).segmentation().type() != "NoSegmentation") {
Sylvester Joosten's avatar
Sylvester Joosten committed
193
          cdim  = converter->cellDimensions(cellID);
Chao Peng's avatar
Chao Peng committed
194
195
196
197
        // get volume dimensions (multiply by two to get fullsize)
        } else {
          // cdim = converter->findContext(id)->volumePlacement().volume().solid().dimensions();
          // Using bounding box instead of actual solid so the dimensions are always in dim_x, dim_y, dim_z
Sylvester Joosten's avatar
Sylvester Joosten committed
198
          cdim = converter->findContext(cellID)->volumePlacement().volume().boundingBox().dimensions();
Chao Peng's avatar
Chao Peng committed
199
200
201
          std::transform(cdim.begin(), cdim.end(), cdim.begin(),
               std::bind(std::multiplies<double>(), std::placeholders::_1, 2));
        }
Sylvester Joosten's avatar
Sylvester Joosten committed
202
203
204
205
206
207
208
209
210
        double dim[3] = {0., 0., 0.};
        for (size_t i = 0; i < cdim.size() && i < 3; ++i) {
          dim[i] = cdim[i] / m_lUnit;
        }
        // debug() << std::bitset<64>(id) << "\n"
        //         <<
        //         m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().volIDs().str()
        //         << endmsg;
        hits.push_back({
211
212
213
214
215
216
217
218
219
220
            rh.cellID(),  // cellID
            rh.ID(),      // ID
            lid,          // layer
            sid,          // sector
            0,            // @TODO: hit type
            energy,       // energy
            0,            // @TODO: energy error
            time,         // time
            {gpos.x() / m_lUnit, gpos.y() / m_lUnit, gpos.z() / m_lUnit}, // global pos
            {pos.x() / m_lUnit, pos.y() / m_lUnit, pos.z() / m_lUnit},    // local pos
Sylvester Joosten's avatar
Sylvester Joosten committed
221
            {dim[0], dim[1], dim[2]}
Sylvester Joosten's avatar
Sylvester Joosten committed
222
223
224
225
        });
      }

      return StatusCode::SUCCESS;
226
227
    }

Sylvester Joosten's avatar
Sylvester Joosten committed
228
  }; // class CalorimeterHitReco
229

Sylvester Joosten's avatar
Sylvester Joosten committed
230
  DECLARE_COMPONENT(CalorimeterHitReco)
231
232

} // namespace Jug::Reco