ATLAS Offline Software
DefectsEmulatorAlg.icc
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 #include "Identifier/Identifier.h"
5 #include "AtlasDetDescr/AtlasDetectorID.h"
6 #include "EventContainers/IdentifiableContainerBase.h"
7 #include "InDetRawData/InDetRawDataCollection.h"
8 #include "InDetRawData/SCT3_RawData.h"
9 #include "InDetRawData/Pixel1RawData.h"
10 #include "InDetReadoutGeometry/SiDetectorElement.h"
11 
12 #include "AthenaKernel/RNGWrapper.h"
13 #include "CLHEP/Random/RandPoisson.h"
14 #include "CLHEP/Random/RandFlat.h"
15 
16 #include "StoreGate/WriteHandle.h"
17 
18 #include <vector>
19 #include <memory>
20 #include <mutex>
21 
22 namespace InDet{
23 
24  template <class T_RDO_Container>
25  bool DefectsEmulatorAlg<T_RDO_Container>::setModuleData(const ActsDetectorElement &acts_detector_element,
26  ModuleIdentifierMatchUtil::ModuleData_t &module_data) const {
27  // It is assumed that a module is unambiguously identified by the detector type and the id hash,
28  // and all modules with the same detector type are relevant for the given RDO container.
29  if (acts_detector_element.detectorType() == DETECTOR_TYPE) {
30  const auto*si_detector_element = dynamic_cast<const InDetDD::SiDetectorElement*>(acts_detector_element.upstreamDetectorElement());
31  if(si_detector_element != nullptr) {
32  const T_ModuleDesign &moduleDesign = dynamic_cast<const T_ModuleDesign &>(si_detector_element->design());
33  ModuleIdentifierMatchUtil::setModuleData(*m_idHelper,
34  acts_detector_element.identify(),
35  moduleDesign,
36  module_data);
37  return true;
38  }
39  }
40  return false;
41  }
42 
43  template <class T_RDO_Container>
44  StatusCode DefectsEmulatorAlg<T_RDO_Container>::initialize(){
45  ATH_CHECK( m_emulatedDefects.initialize() );
46  ATH_CHECK( m_rdoOutContainerKey.initialize() );
47  ATH_CHECK( m_origRdoContainerKey.initialize() );
48  ATH_CHECK( detStore()->retrieve(m_idHelper, m_idHelperName.value()) );
49 
50  return DefectsEmulatorBase::initializeBase(m_idHelper->wafer_hash_max());
51  }
52 
53  template <class T_RDO_Container>
54  StatusCode DefectsEmulatorAlg<T_RDO_Container>::execute(const EventContext& ctx) const {
55  SG::ReadCondHandle<T_DefectsData> emulatedDefects(m_emulatedDefects,ctx);
56  ATH_CHECK(emulatedDefects.isValid());
57  SG::ReadHandle<T_RDO_Container> origRdoContainer(m_origRdoContainerKey, ctx);
58  ATH_CHECK(origRdoContainer.isValid());
59  SG::WriteHandle<T_RDO_Container> rdoOutContainer(m_rdoOutContainerKey, ctx);
60  ATH_CHECK( rdoOutContainer.record (std::make_unique<T_RDO_Container>(origRdoContainer->size(), EventContainers::Mode::OfflineFast)) );
61 
62  CLHEP::HepRandomEngine *rndmEngine{};
63  if (!m_noiseProbability.value().empty()){
64  m_rndmSvc->getEngine(this, m_rngName);
65  ATHRNG::RNGWrapper* rngWrapper = m_rndmSvc->getEngine(this, m_rngName);
66  rngWrapper->setSeed( m_rngName, ctx );
67  rndmEngine = rngWrapper->getEngine(ctx);
68  }
69 
70  unsigned int n_rejected=0u;
71  unsigned int n_split_rdos=0u;
72  unsigned int n_new=0;
73  T_ID_Adapter id_helper(m_idHelper);
74  unsigned int n_added_noise_hits=0u;
75  unsigned int n_defect_noise_hits=0u;
76 
77  std::vector<typename T_ModuleHelper::KEY_TYPE> sorted_keys;
78  for(const InDetRawDataCollection<T_RDORawData>* collection : *origRdoContainer ) {
79  const IdentifierHash idHash = collection->identifyHash();
80  T_ModuleHelper module_helper( emulatedDefects->getDetectorElement(idHash).design() );
81  std::unique_ptr<InDetRawDataCollection<T_RDORawData> >
82  clone = std::make_unique<InDetRawDataCollection<T_RDORawData> >(idHash);
83  clone->setIdentifier(collection->identify());
84  unsigned int rejected_per_mod=0u;
85  unsigned int noise_per_mod=0u;
86  if (!emulatedDefects->isModuleDefect(idHash)) {
87  unsigned int n_cells = module_helper.nCells();
88  unsigned int n_noise_hits=rndmEngine && m_noiseParamIdx.at(idHash)<m_noiseProbability.size()
89  ? CLHEP::RandPoisson::shoot(rndmEngine,
90  n_cells * m_noiseProbability.value()[ m_noiseParamIdx[idHash] ] )
91  : 0u;
92  clone->reserve( collection->size() + n_noise_hits);
93  if (!module_helper) {
94  ATH_MSG_ERROR( "Not module design for " << idHash);
95  return StatusCode::FAILURE;
96  }
97  TH2 *h2_rejected_hits=nullptr;
98  TH2 *h2_noise_hits=nullptr;
99  TH1 *h1_noise_shape=nullptr;
100  if (m_histogrammingEnabled) {
101  std::lock_guard<std::mutex> lock(m_histMutex);
102  std::tuple<TH2 *, TH2 *, TH1 *> hists = findHist(module_helper.nSensorRows(), module_helper.nSensorColumns());
103  h2_rejected_hits=std::get<0>(hists);
104  h2_noise_hits=std::get<1>(hists);
105  h1_noise_shape=std::get<2>(hists);
106  }
107 
108  sorted_keys.clear();
109  sorted_keys.reserve( collection->size()+n_noise_hits);
110  for(const auto *const rdo : *collection) {
111  const Identifier rdoID = rdo->identify();
112 
113  auto row_idx = id_helper.row_index(rdoID);
114  auto col_idx = id_helper.col_index(rdoID);
115 
116  unsigned int n_new_hits = id_helper.cloneOrRejectHit( module_helper, *(emulatedDefects.cptr()), idHash, row_idx, col_idx, *rdo, *clone);
117  if (n_new_hits>0) {
118  n_split_rdos += (n_new_hits-1);
119  for (unsigned int hit_i=n_new_hits; hit_i>0; --hit_i) {
120  auto new_hit = clone->at(clone->size()-hit_i);
121 
122  const Identifier rdoID = new_hit->identify();
123  auto row_idx = id_helper.row_index(rdoID);
124  auto col_idx = id_helper.col_index(rdoID);
125  unsigned int key = module_helper.hardwareCoordinates(row_idx,col_idx);
126  auto iter = std::lower_bound(sorted_keys.begin(), sorted_keys.end(), key);
127  if (iter == sorted_keys.end()) {
128  sorted_keys.push_back(key);
129  }
130  else {
131  if (*iter != key) {
132  sorted_keys.insert(iter, key);
133  }
134  }
135  }
136  }
137  else {
138  if (m_histogrammingEnabled && h2_rejected_hits) {
139  std::lock_guard<std::mutex> lock(m_histMutex);
140  h2_rejected_hits->Fill(col_idx, row_idx);
141  }
142  ++rejected_per_mod;
143  }
144 
145  }
146 
147  // add random noise hits
148  for (unsigned int i=0; i<n_noise_hits; ++i) {
149  for (unsigned int attempt_i=0; attempt_i<10; ++attempt_i) {
150  unsigned int cell_idx=CLHEP::RandFlat::shoot(rndmEngine,n_cells); // %pixels;
151  unsigned int row_idx = cell_idx % module_helper.rows();
152  unsigned int col_idx = cell_idx / module_helper.rows();
153  unsigned int key = module_helper.hardwareCoordinates(row_idx,col_idx);
154  if (emulatedDefects->isDefect(module_helper, idHash, row_idx, col_idx)) {
155  ++n_defect_noise_hits;
156  break;
157  }
158  auto iter = std::lower_bound(sorted_keys.begin(), sorted_keys.end(), key);
159  if (iter == sorted_keys.end()) {
160  sorted_keys.push_back(key);
161  }
162  else {
163  if (*iter != key) {
164  sorted_keys.insert(iter, key);
165  }
166  else {
167  continue;
168  }
169  }
170  float prob = CLHEP::RandFlat::shoot(rndmEngine,1.);
171  assert( m_noiseParamIdx[idHash] < m_noiseShapeCummulative.size() );
172  const std::vector<float> &noise_shape = m_noiseShapeCummulative[ m_noiseParamIdx[idHash] ];
173  unsigned int noise_value=0;
174  for (noise_value=0; noise_value<noise_shape.size(); ++noise_value) {
175  if (noise_shape[noise_value]>prob) {
176  break;
177  }
178  }
179 
180  clone->push_back( id_helper.createNoiseHit( module_helper, collection->identify(), cell_idx, noise_value).release() );
181  if (m_histogrammingEnabled && h2_noise_hits) {
182  unsigned int row_idx = cell_idx % module_helper.rows();
183  unsigned int col_idx = cell_idx / module_helper.rows();
184 
185  std::lock_guard<std::mutex> lock(m_histMutex);
186  h2_noise_hits->Fill(col_idx, row_idx);
187  h1_noise_shape->Fill(noise_value);
188  }
189  ++noise_per_mod;
190  break;
191  }
192  }
193 
194 
195  }
196  else {
197  rejected_per_mod += collection->size();
198  }
199  if (m_histogrammingEnabled) {
200  std::lock_guard<std::mutex> lock(m_histMutex);
201  unsigned int bin_i=m_moduleHist[kRejectedHits]->GetBin( idHash%100+1, idHash/100+1);
202  assert(m_moduleHist[kRejectedHits]->GetNbinsX() == m_moduleHist[kNoiseHits]->GetNbinsX());
203  assert(m_moduleHist[kRejectedHits]->GetNbinsY() == m_moduleHist[kNoiseHits]->GetNbinsY());
204  m_moduleHist[kRejectedHits]->SetBinContent(bin_i, rejected_per_mod );
205  m_moduleHist[kNoiseHits]->SetBinContent(bin_i, noise_per_mod );
206  }
207  n_added_noise_hits += noise_per_mod;
208  n_rejected += rejected_per_mod;
209  n_new += clone->size();
210  rdoOutContainer->addCollection( clone.release(), idHash).ignore();
211  }
212  m_rejectedRDOs += n_rejected;
213  m_totalRDOs += n_new;
214  m_totalNoise += n_added_noise_hits;
215  m_splitRDOs += n_split_rdos;
216  ATH_MSG_DEBUG("rejected " << m_rejectedRDOs << ", copied " << m_totalRDOs << " RDOs"
217  << " noise " << n_added_noise_hits << " (+" << n_defect_noise_hits << ")"
218  << " split hits: " << n_split_rdos);
219 
220  return StatusCode::SUCCESS;
221  }
222 
223 }