2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
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"
12 #include "AthenaKernel/RNGWrapper.h"
13 #include "CLHEP/Random/RandPoisson.h"
14 #include "CLHEP/Random/RandFlat.h"
16 #include "StoreGate/WriteHandle.h"
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(),
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()) );
50 return DefectsEmulatorBase::initializeBase(m_idHelper->wafer_hash_max());
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)) );
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);
70 unsigned int n_rejected=0u;
71 unsigned int n_split_rdos=0u;
73 T_ID_Adapter id_helper(m_idHelper);
74 unsigned int n_added_noise_hits=0u;
75 unsigned int n_defect_noise_hits=0u;
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] ] )
92 clone->reserve( collection->size() + n_noise_hits);
94 ATH_MSG_ERROR( "Not module design for " << idHash);
95 return StatusCode::FAILURE;
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);
109 sorted_keys.reserve( collection->size()+n_noise_hits);
110 for(const auto *const rdo : *collection) {
111 const Identifier rdoID = rdo->identify();
113 auto row_idx = id_helper.row_index(rdoID);
114 auto col_idx = id_helper.col_index(rdoID);
116 unsigned int n_new_hits = id_helper.cloneOrRejectHit( module_helper, *(emulatedDefects.cptr()), idHash, row_idx, col_idx, *rdo, *clone);
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);
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);
132 sorted_keys.insert(iter, key);
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);
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;
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);
164 sorted_keys.insert(iter, key);
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) {
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();
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);
197 rejected_per_mod += collection->size();
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 );
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();
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);
220 return StatusCode::SUCCESS;