CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
HGCalSoARecHitsProducer.cc
Go to the documentation of this file.
16 
19 
21 
23  public:
25  : detector_(config.getParameter<std::string>("detector")),
27  isNose_(detector_ == "HFNose"),
28  maxNumberOfThickIndices_(config.getParameter<unsigned>("maxNumberOfThickIndices")),
29  fcPerEle_(config.getParameter<double>("fcPerEle")),
30  ecut_(config.getParameter<double>("ecut")),
31  fcPerMip_(config.getParameter<std::vector<double>>("fcPerMip")),
32  nonAgedNoises_(config.getParameter<std::vector<double>>("noises")),
33  dEdXweights_(config.getParameter<std::vector<double>>("dEdXweights")),
34  thicknessCorrection_(config.getParameter<std::vector<double>>("thicknessCorrection")),
36  hits_token_(consumes<HGCRecHitCollection>(config.getParameter<edm::InputTag>("recHits"))),
38 
39  ~HGCalSoARecHitsProducer() override = default;
40 
41  void produce(device::Event& iEvent, device::EventSetup const& iSetup) override {
43 
47 
48  hits_h = iEvent.getHandle(hits_token_);
49  auto const& hits = *(hits_h.product());
51 
52  // Count effective hits above threshold
53  uint32_t index = 0;
54  for (unsigned int i = 0; i < hits.size(); ++i) {
55  const HGCRecHit& hgrh = hits[i];
56  DetId detid = hgrh.detid();
57  unsigned int layerOnSide = (rhtools_.getLayerWithOffset(detid) - 1);
58 
59  // set sigmaNoise default value 1 to use kappa value directly in case of
60  // sensor-independent thresholds
61  int thickness_index = rhtools_.getSiThickIndex(detid);
62  if (thickness_index == -1) {
63  thickness_index = maxNumberOfThickIndices_;
64  }
65  double storedThreshold = thresholds_[layerOnSide][thickness_index];
66  if (hgrh.energy() < storedThreshold)
67  continue; // this sets the ZS threshold at ecut times the sigma noise
68  index++;
69  }
70 
71  // Allocate Host SoA will contain one entry for each RecHit above threshold
73  auto cellsView = cells.view();
74 
75  // loop over all hits and create the Hexel structure, skip energies below ecut
76  // for each layer and wafer calculate the thresholds (sigmaNoise and energy)
77  // once
78  index = 0;
79  for (unsigned int i = 0; i < hits.size(); ++i) {
80  const HGCRecHit& hgrh = hits[i];
81  DetId detid = hgrh.detid();
82  unsigned int layerOnSide = (rhtools_.getLayerWithOffset(detid) - 1);
83 
84  // set sigmaNoise default value 1 to use kappa value directly in case of
85  // sensor-independent thresholds
86  float sigmaNoise = 1.f;
87  int thickness_index = rhtools_.getSiThickIndex(detid);
88  if (thickness_index == -1) {
89  thickness_index = maxNumberOfThickIndices_;
90  }
91  double storedThreshold = thresholds_[layerOnSide][thickness_index];
92  if (detid.det() == DetId::HGCalHSi || detid.subdetId() == HGCHEF) {
93  storedThreshold = thresholds_.at(layerOnSide).at(thickness_index + deltasi_index_regemfac_);
94  }
95  sigmaNoise = v_sigmaNoise_.at(layerOnSide).at(thickness_index);
96 
97  if (hgrh.energy() < storedThreshold)
98  continue; // this sets the ZS threshold at ecut times the sigma noise
99  // for the sensor
100 
102  int offset = ((rhtools_.zside(detid) + 1) >> 1) * maxlayer_;
103  int layer = layerOnSide + offset;
104  auto entryInSoA = cellsView[index];
105  if (detector_ == "BH") {
106  entryInSoA.dim1() = position.eta();
107  entryInSoA.dim2() = position.phi();
108  } // else, isSilicon == true and eta phi values will not be used
109  else {
110  entryInSoA.dim1() = position.x();
111  entryInSoA.dim2() = position.y();
112  }
113  entryInSoA.dim3() = position.z();
114  entryInSoA.weight() = hgrh.energy();
115  entryInSoA.sigmaNoise() = sigmaNoise;
116  entryInSoA.layer() = layer;
117  entryInSoA.recHitIndex() = i;
118  entryInSoA.detid() = detid.rawId();
119  entryInSoA.time() = hgrh.time();
120  entryInSoA.timeError() = hgrh.timeError();
121  index++;
122  }
123 #if 0
124  std::cout << "Size: " << cells->metadata().size() << " count cells: " << index
125  << " i.e. " << cells->metadata().size() << std::endl;
126 #endif
127 
128  if constexpr (!std::is_same_v<Device, alpaka_common::DevHost>) {
129  // Trigger copy async to GPU
130  //std::cout << "GPU" << std::endl;
131  HGCalSoARecHitsDeviceCollection deviceProduct{cells->metadata().size(), iEvent.queue()};
132  alpaka::memcpy(iEvent.queue(), deviceProduct.buffer(), cells.const_buffer());
133  iEvent.emplace(deviceToken_, std::move(deviceProduct));
134  } else {
135  //std::cout << "CPU" << std::endl;
136  iEvent.emplace(deviceToken_, std::move(cells));
137  }
138  }
139 
142  desc.add<std::string>("detector", "EE")->setComment("options EE, FH, BH, HFNose; other value defaults to EE");
143  desc.add<edm::InputTag>("recHits", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
144  desc.add<unsigned int>("maxNumberOfThickIndices", 6);
145  desc.add<double>("fcPerEle", 0.00016020506);
146  desc.add<std::vector<double>>("fcPerMip");
147  desc.add<std::vector<double>>("thicknessCorrection");
148  desc.add<std::vector<double>>("noises");
149  desc.add<std::vector<double>>("dEdXweights");
150  desc.add<double>("ecut", 3.);
151  descriptions.addWithDefaultLabel(desc);
152  }
153 
154  private:
157  bool isNose_;
159  unsigned int maxlayer_;
162  double fcPerEle_;
163  double ecut_;
164  std::vector<double> fcPerMip_;
165  std::vector<double> nonAgedNoises_;
166  std::vector<double> dEdXweights_;
167  std::vector<double> thicknessCorrection_;
168  std::vector<std::vector<double>> thresholds_;
169  std::vector<std::vector<double>> v_sigmaNoise_;
170 
175 
177  // To support the TDR geometry and also the post-TDR one (v9 onwards), we
178  // need to change the logic of the vectors containing signal to noise and
179  // thresholds. The first 3 indices will keep on addressing the different
180  // thicknesses of the Silicon detectors in CE_E , the next 3 indices will address
181  // the thicknesses of the Silicon detectors in CE_H, while the last one, number 6 (the
182  // seventh) will address the Scintillators. This change will support both
183  // geometries at the same time.
184 
185  if (initialized_)
186  return; // only need to calculate thresholds once
187 
188  initialized_ = true;
189 
190  std::vector<double> dummy;
191 
192  dummy.resize(maxNumberOfThickIndices_ + !isNose_, 0); // +1 to accomodate for the Scintillators
193  thresholds_.resize(maxlayer_, dummy);
194  v_sigmaNoise_.resize(maxlayer_, dummy);
195 
196  for (unsigned ilayer = 1; ilayer <= maxlayer_; ++ilayer) {
197  for (unsigned ithick = 0; ithick < maxNumberOfThickIndices_; ++ithick) {
198  float sigmaNoise = 0.001f * fcPerEle_ * nonAgedNoises_[ithick] * dEdXweights_[ilayer] /
199  (fcPerMip_[ithick] * thicknessCorrection_[ithick]);
200  thresholds_[ilayer - 1][ithick] = sigmaNoise * ecut_;
201  v_sigmaNoise_[ilayer - 1][ithick] = sigmaNoise;
202 #if 0
203  std::cout << "ilayer: " << ilayer << " nonAgedNoises: " << nonAgedNoises_[ithick]
204  << " fcPerEle: " << fcPerEle_ << " fcPerMip: " << fcPerMip_[ithick]
205  << " noiseMip: " << fcPerEle_ * nonAgedNoises_[ithick] / fcPerMip_[ithick]
206  << " sigmaNoise: " << sigmaNoise << "\n";
207 #endif
208  }
209  }
210  }
211  };
212 
213 } // namespace ALPAKA_ACCELERATOR_NAMESPACE
214 
216 DEFINE_FWK_ALPAKA_MODULE(HGCalSoARecHitsProducer);
device::EDPutToken< HGCalSoARecHitsDeviceCollection > const deviceToken_
void addWithDefaultLabel(ParameterSetDescription const &psetDescription)
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
std::vector< HGCRecHit > HGCRecHitCollection
constexpr const DetId & detid() const
Definition: CaloRecHit.h:33
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
T const * product() const
Definition: Handle.h:70
edm::ESHandle< T > getHandle(edm::ESGetToken< T, R > const &iToken) const
Definition: EventSetup.h:49
Definition: config.py:1
int zside(const DetId &id) const
Definition: RecHitTools.cc:174
constexpr float energy() const
Definition: CaloRecHit.h:29
PortableCollection< HGCalSoARecHits > HGCalSoARecHitsDeviceCollection
int iEvent
Definition: GenABIO.cc:224
GlobalPoint getPosition(const DetId &id) const
Definition: RecHitTools.cc:140
void produce(device::Event &iEvent, device::EventSetup const &iSetup) override
Definition: DetId.h:17
void setGeometry(CaloGeometry const &)
Definition: RecHitTools.cc:79
constexpr float time() const
Definition: CaloRecHit.h:31
HLT enums.
edm::ESGetToken< CaloGeometry, CaloGeometryRecord > caloGeomToken_
static int position[264][3]
Definition: ReadPGInfo.cc:289
float timeError() const
Definition: HGCRecHit.cc:79
auto produces(std::string instanceName) noexcept
declare what type of product will make and with which optional label
#define DEFINE_FWK_ALPAKA_MODULE(name)
Definition: MakerMacros.h:16
int getSiThickIndex(const DetId &) const
Definition: RecHitTools.cc:216
def move(src, dest)
Definition: eostools.py:511
unsigned int lastLayer(bool nose=false) const
Definition: RecHitTools.h:80
unsigned int getLayerWithOffset(const DetId &) const
Definition: RecHitTools.cc:381