ATLAS Offline Software
Loading...
Searching...
No Matches
RpcDigitizationTool.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
6//
7// RpcDigitizationTool
8// ------------
9// Authors:
10// Andrea Di Simone <Andrea.Di.Simone@cern.ch>
11// Gabriele Chiodini <gabriele.chiodini@le.infn.it>
12// Stefania Spagnolo <stefania.spagnolo@le.infn.it>
14
15#include "RPC_Digitization/RpcDigitizationTool.h"
16
17// Inputs
18#include "GaudiKernel/SystemOfUnits.h"
19#include "GaudiKernel/PhysicalConstants.h"
20#include "GeoPrimitives/GeoPrimitivesToStringConverter.h"
21#include "GeoModelHelpers/TransformToStringConverter.h"
22#include "MuonSimEvent/RPCSimHit.h"
23#include "MuonSimEvent/RPCSimHitCollection.h"
24
25// Geometry
26#include "MuonIdHelpers/RpcIdHelper.h"
27#include "MuonReadoutGeometry/MuonDetectorManager.h"
28#include "MuonReadoutGeometry/RpcReadoutElement.h"
29#include "MuonSimEvent/RpcHitIdHelper.h"
30
31// run n. from geometry DB
32#include "GeoModelInterfaces/IGeoModelSvc.h"
33#include "RDBAccessSvc/IRDBAccessSvc.h"
34#include "RDBAccessSvc/IRDBRecord.h"
35#include "RDBAccessSvc/IRDBRecordset.h"
36
37// Truth
38#include "AtlasHepMC/GenParticle.h"
39#include "TruthUtils/HepMCHelpers.h"
40#include "GeneratorObjects/HepMcParticleLink.h"
41#include "GeoModelKernel/throwExcept.h"
42// Random Numbers
43#include "AthenaKernel/RNGWrapper.h"
44#include "CLHEP/Random/RandExponential.h"
45#include "CLHEP/Random/RandFlat.h"
46#include "CLHEP/Random/RandGaussZiggurat.h"
47
48// Core includes
49#include <TString.h> // for Form
50
51#include <atomic>
52#include <fstream>
53#include <iostream>
54#include <sstream>
55#include <utility>
56
57#include "EventInfoMgt/ITagInfoMgr.h"
58#include "PathResolver/PathResolver.h"
59
60// 12 charge points, 15 BetaGamma points, 180 efficiency points for fcp search
61namespace {
62 constexpr int N_Charge = 12;
63 constexpr int N_Velocity = 15;
64 constexpr std::array<double, N_Charge> Charge{0.1, 0.2, 0.3, 0.33, 0.4, 0.5, 0.6, 0.66, 0.7, 0.8, 0.9, 1.0};
65 constexpr std::array<double, N_Velocity> Velocity{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 2.0, 3.0, 10.0, 100.0, 1000.0};
66 constexpr double Eff_garfield[N_Charge][N_Velocity] = {
67 {0.8648, 0.3476, 0.1407, 0.0618, 0.0368, 0.0234, 0.0150, 0.0120, 0.0096, 0.0079, 0.0038, 0.0041, 0.0035, 0.0049, 0.0054},
68 {0.9999, 0.9238, 0.6716, 0.4579, 0.3115, 0.2238, 0.1727, 0.1365, 0.1098, 0.0968, 0.0493, 0.0451, 0.0528, 0.0694, 0.0708},
69 {1.0000, 0.9978, 0.9517, 0.8226, 0.6750, 0.5611, 0.4674, 0.3913, 0.3458, 0.3086, 0.1818, 0.1677, 0.1805, 0.2307, 0.2421},
70 {1.0000, 0.9994, 0.9758, 0.8918, 0.7670, 0.6537, 0.5533, 0.4856, 0.4192, 0.3852, 0.2333, 0.2186, 0.2479, 0.2957, 0.2996},
71 {1.0000, 1.0000, 0.9972, 0.9699, 0.9022, 0.8200, 0.7417, 0.6660, 0.6094, 0.5622, 0.3846, 0.3617, 0.3847, 0.4578, 0.4583},
72 {1.0000, 1.0000, 0.9998, 0.9956, 0.9754, 0.9479, 0.9031, 0.8604, 0.8126, 0.7716, 0.5827, 0.5545, 0.5865, 0.6834, 0.6706},
73 {1.0000, 1.0000, 1.0000, 0.9997, 0.9968, 0.9876, 0.9689, 0.9464, 0.9221, 0.8967, 0.7634, 0.7385, 0.7615, 0.8250, 0.8309},
74 {1.0000, 1.0000, 1.0000, 1.0000, 0.9995, 0.9952, 0.9866, 0.9765, 0.9552, 0.9427, 0.8373, 0.8127, 0.8412, 0.8899, 0.8891},
75 {1.0000, 1.0000, 1.0000, 1.0000, 0.9995, 0.9981, 0.9918, 0.9803, 0.9754, 0.9602, 0.8730, 0.8564, 0.8746, 0.9178, 0.9261},
76 {1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 0.9993, 0.9990, 0.9951, 0.9935, 0.9886, 0.9419, 0.9277, 0.9422, 0.9686, 0.9700},
77 {1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 0.9998, 0.9996, 0.9980, 0.9966, 0.9786, 0.9718, 0.9748, 0.9875, 0.9882},
78 {1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 0.9998, 1.0000, 0.9991, 0.9988, 0.9913, 0.9872, 0.9917, 0.9970, 0.9964}};
79 bool
80 validIndex(int idx, int arraySize){
81 return (idx>=0) and (idx<arraySize);
82 }
83} // namespace
84
85using namespace MuonGM;
86namespace {
87 constexpr double SIG_VEL = 4.8;
88}
89
90
91RpcDigitizationTool::RpcDigitizationTool(const std::string& type, const std::string& name, const IInterface* pIID) :
92 PileUpToolBase(type, name, pIID) {}
93
94// member function implementation
95//--------------------------------------------
96StatusCode RpcDigitizationTool::initialize() {
97 ATH_MSG_DEBUG("RpcDigitizationTool:: in initialize()");
98 ATH_MSG_DEBUG("Configuration RpcDigitizationTool ");
99
100 ATH_MSG_DEBUG("InputObjectName " << m_inputHitCollectionName);
101 ATH_MSG_DEBUG("OutputObjectName " << m_outputDigitCollectionKey.key());
102 ATH_MSG_DEBUG("OutputSDOName " << m_outputSDO_CollectionKey.key());
103 ATH_MSG_DEBUG("WindowLowerOffset " << m_timeWindowLowerOffset);
104 ATH_MSG_DEBUG("WindowUpperOffset " << m_timeWindowUpperOffset);
105 ATH_MSG_DEBUG("DeadTime " << m_deadTime);
106 ATH_MSG_DEBUG("RndmSvc " << m_rndmSvc);
107 ATH_MSG_DEBUG("PatchForRpcTime " << m_patch_for_rpc_time);
108 ATH_MSG_DEBUG("RpcTimeShift " << m_rpc_time_shift);
109 ATH_MSG_DEBUG("RPC_TimeSchema " << m_RPC_TimeSchema);
110 ATH_MSG_DEBUG("RPCSDOareRPCDigits " << m_sdoAreOnlyDigits);
111
112 ATH_MSG_DEBUG("IgnoreRunDependentConfig " << m_ignoreRunDepConfig);
113 ATH_MSG_DEBUG("turnON_efficiency " << m_turnON_efficiency);
114 ATH_MSG_DEBUG("Efficiency_fromCOOL " << m_Efficiency_fromCOOL);
115 ATH_MSG_DEBUG("Efficiency_BIS78_fromCOOL" << m_Efficiency_BIS78_fromCOOL);
116 ATH_MSG_DEBUG("turnON_clustersize " << m_turnON_clustersize);
117 ATH_MSG_DEBUG("ClusterSize_fromCOOL " << m_ClusterSize_fromCOOL);
118 ATH_MSG_DEBUG("ClusterSize_BIS78_fromCOOL" << m_ClusterSize_BIS78_fromCOOL);
119 ATH_MSG_DEBUG("FirstClusterSizeInTail " << m_FirstClusterSizeInTail);
120 ATH_MSG_DEBUG("ClusterSize1_2uncorr " << m_ClusterSize1_2uncorr);
121 ATH_MSG_DEBUG("BOG_BOF_DoubletR2_OFF " << m_BOG_BOF_DoubletR2_OFF);
122 ATH_MSG_DEBUG("CutMaxClusterSize " << m_CutMaxClusterSize);
123 ATH_MSG_DEBUG("CutProjectedTracks " << m_CutProjectedTracks);
124 ATH_MSG_DEBUG("ValidationSetup " << m_validationSetup);
125 ATH_MSG_DEBUG("IncludePileUpTruth " << m_includePileUpTruth);
126 ATH_MSG_DEBUG("VetoPileUpTruthLinks " << m_vetoPileUpTruthLinks);
127
128 ATH_CHECK(m_detMgrKey.initialize());
129 if (m_onlyUseContainerName) { ATH_CHECK(m_mergeSvc.retrieve()); }
130 ATH_CHECK(detStore()->retrieve(m_idHelper));
131 // check the identifiers
132
133 ATH_MSG_INFO("Max Number of RPC Gas Gaps for these Identifiers = " << m_idHelper->gasGapMax());
134
135 // check the input object name
136 if (m_hitsContainerKey.key().empty()) {
137 ATH_MSG_FATAL("Property InputObjectName not set !");
138 return StatusCode::FAILURE;
139 }
140 if (m_onlyUseContainerName) m_inputHitCollectionName = m_hitsContainerKey.key();
141 ATH_MSG_DEBUG("Input objects in container : '" << m_inputHitCollectionName << "'");
142
143 // Initialize ReadHandleKey
144 ATH_CHECK(m_hitsContainerKey.initialize());
145
146 // initialize the output WriteHandleKeys
147 ATH_CHECK(m_outputDigitCollectionKey.initialize());
148 ATH_CHECK(m_outputSDO_CollectionKey.initialize());
149 ATH_CHECK(m_simHitValidKey.initialize(m_validationSetup));
150 ATH_MSG_DEBUG("Output digits: '" << m_outputDigitCollectionKey.key() << "'");
151
152 // set the configuration based on run1/run2
153 ATH_CHECK(initializeRunDependentParameters());
154
155 ATH_MSG_DEBUG("Ready to read parameters for cluster simulation from file");
156
157 ATH_CHECK(m_rndmSvc.retrieve());
158
159 // fill the taginfo information
160 ATH_CHECK(fillTagInfo());
161
162 ATH_CHECK(m_readKey.initialize(m_RPCInfoFromDb));
163
165 // m_turnON_clustersize=false;
166 m_BOF_id = m_idHelper->stationNameIndex("BOF");
167 m_BOG_id = m_idHelper->stationNameIndex("BOG");
168 m_BOS_id = m_idHelper->stationNameIndex("BOS");
169 m_BIL_id = m_idHelper->stationNameIndex("BIL");
170 m_BIS_id = m_idHelper->stationNameIndex("BIS");
171 m_muonHelper = RpcHitIdHelper::GetHelper(m_idHelper->gasGapMax());
172
173 return StatusCode::SUCCESS;
174}
175
176StatusCode RpcDigitizationTool::initializeRunDependentParameters() {
177 // TODO This should all be in a conditions Alg
178 // Retrieve geometry config information from the database (RUN1, RUN2, etc...)
179 SmartIF<IGeoModelSvc> geoModel{Gaudi::svcLocator()->service("GeoModelSvc")};
180 if ( !geoModel ) {
181 ATH_MSG_ERROR("Could not locate GeoModelSvc");
182 return StatusCode::FAILURE;
183 }
184
185 // check the DetDescr version
186 std::string atlasVersion = geoModel->atlasVersion();
187
188 SmartIF<IRDBAccessSvc> rdbAccess{Gaudi::svcLocator()->service("RDBAccessSvc")};
189 if ( !rdbAccess ) {
190 ATH_MSG_ERROR("Could not locate RDBAccessSvc");
191 return StatusCode::FAILURE;
192 }
193
194 enum DataPeriod {Unknown, Run1, Run2, Run3, Run4 };
195 DataPeriod run = Unknown;
196
197 std::string configVal = "";
198
199 IRDBRecordset_ptr atlasCommonRec = rdbAccess->getRecordsetPtr("AtlasCommon", atlasVersion, "ATLAS");
200 if (atlasCommonRec->size() == 0) {
201 run = Run1;
202 } else {
203 configVal = (*atlasCommonRec)[0]->getString("CONFIG");
204 ATH_MSG_INFO("From DD Database, Configuration is " << configVal);
205 if (configVal == "RUN1") {
206 run = Run1;
207 } else if (configVal == "RUN2") {
208 run = Run2;
209 } else if (configVal == "RUN3") {
210 run = Run3;
211 } else if (configVal == "RUN4") {
212 run = Run4;
213 }
214 if (run == DataPeriod::Unknown) {
215 ATH_MSG_FATAL("Unexpected value for geometry config read from the database: " << configVal);
216 return StatusCode::FAILURE;
217 }
218 }
219 if (run == Run3 && m_idHelper->gasGapMax() < 3)
220 ATH_MSG_WARNING("Run3, configVal = " << configVal << " and GasGapMax =" << m_idHelper->gasGapMax());
221
222 if (run == Run1)
223 ATH_MSG_INFO("From Geometry DB: MuonSpectrometer configuration is: RUN1 or MuonGeometry = R.06");
224 else if (run == Run2)
225 ATH_MSG_INFO("From Geometry DB: MuonSpectrometer configuration is: RUN2 or MuonGeometry = R.07");
226 else if (run == Run3)
227 ATH_MSG_INFO("From Geometry DB: MuonSpectrometer configuration is: RUN3 or MuonGeometry = R.09");
228 else if (run == Run4)
229 ATH_MSG_INFO("From Geometry DB: MuonSpectrometer configuration is: RUN4 or MuonGeometry = R.10");
230
231 if (m_ignoreRunDepConfig == false) {
232 m_BOG_BOF_DoubletR2_OFF = false;
233 m_Efficiency_fromCOOL = false;
234 m_ClusterSize_fromCOOL = false;
235 m_RPCInfoFromDb = false;
236 m_kill_deadstrips = false;
237 if (run == Run1) {
238 // m_BOG_BOF_DoubletR2_OFF = true
239 // m_Efficiency_fromCOOL = true
240 // m_ClusterSize_fromCOOL = true
241 m_BOG_BOF_DoubletR2_OFF = true;
242 if (configVal == "RUN1") { // MC12 setup
243 m_Efficiency_fromCOOL = true;
244 m_ClusterSize_fromCOOL = true;
245 m_RPCInfoFromDb = true;
246 m_kill_deadstrips = true;
247 m_CutProjectedTracks = 50;
248 }
249 } else {
250 // m_BOG_BOF_DoubletR2_OFF = false # do not turn off at digitization the hits in the dbR=2 chambers in the feet
251 // m_Efficiency_fromCOOL = false # use common average values in python conf.
252 // m_ClusterSize_fromCOOL = false # use common average values in python conf.
253 m_BOG_BOF_DoubletR2_OFF = false;
254 if (run == Run2) { // MC15c setup
255 m_Efficiency_fromCOOL = true;
256 m_ClusterSize_fromCOOL = true;
257 m_RPCInfoFromDb = true;
258 m_kill_deadstrips = false;
259 m_CutProjectedTracks = 100;
260 } else {
261 ATH_MSG_INFO("Run3/4: configuration parameter not from COOL");
262 m_Efficiency_fromCOOL = false;
263 m_ClusterSize_fromCOOL = false;
264 m_RPCInfoFromDb = false;
265 m_kill_deadstrips = false;
266 }
267 }
268 ATH_MSG_INFO("RPC Run1/2/3-dependent configuration is enforced");
269 } else {
270 ATH_MSG_WARNING("Run1/2/3-dependent configuration is bypassed; be careful with option settings");
271 }
272
273 ATH_MSG_DEBUG("......RPC Efficiency_fromCOOL " << m_Efficiency_fromCOOL);
274 ATH_MSG_DEBUG("......RPC ClusterSize_fromCOOL " << m_ClusterSize_fromCOOL);
275 ATH_MSG_DEBUG("......RPC BOG_BOF_DoubletR2_OFF " << m_BOG_BOF_DoubletR2_OFF);
276 ATH_MSG_DEBUG("......RPC RPCInfoFromDb " << m_RPCInfoFromDb);
277 ATH_MSG_DEBUG("......RPC KillDeadStrips " << m_kill_deadstrips);
278 ATH_MSG_DEBUG("......RPC CutProjectedTracks " << m_CutProjectedTracks);
279
280
281 return StatusCode::SUCCESS;
282}
283
284template <class CondType>
285StatusCode RpcDigitizationTool::retrieveCondData(const EventContext& ctx,
286 const SG::ReadCondHandleKey<CondType>& key,
287 const CondType* & condPtr) const {
288
289 if (key.empty()) {
290 ATH_MSG_DEBUG("No key has been configured for object "<<typeid(CondType).name()<<". Clear pointer");
291 condPtr = nullptr;
292 return StatusCode::SUCCESS;
293 }
294 SG::ReadCondHandle<CondType> readHandle{key, ctx};
295 if (!readHandle.isValid()){
296 ATH_MSG_FATAL("Failed to load conditions object "<<key.fullKey()<<".");
297 return StatusCode::FAILURE;
298 }
299 condPtr = readHandle.cptr();
300 return StatusCode::SUCCESS;
301
302}
303//--------------------------------------------
304StatusCode RpcDigitizationTool::prepareEvent(const EventContext& /*ctx*/, unsigned int) {
305 ATH_MSG_DEBUG("RpcDigitizationTool::in prepareEvent()");
306
307 // John's Hacks START
308 m_RPCHitCollList.clear();
309 m_thpcRPC = std::make_unique<TimedHitCollection<RPCSimHit>>();
310 // John's Hacks END
311
312 return StatusCode::SUCCESS;
313}
314
315//--------------------------------------------
316StatusCode RpcDigitizationTool::processBunchXing(int bunchXing, SubEventIterator bSubEvents, SubEventIterator eSubEvents) {
317 ATH_MSG_DEBUG("RpcDigitizationTool::in processBunchXing()");
318
319 typedef PileUpMergeSvc::TimedList<RPCSimHitCollection>::type TimedHitCollList;
320 TimedHitCollList hitCollList;
321
322 if (!(m_mergeSvc->retrieveSubSetEvtData(m_inputHitCollectionName, hitCollList, bunchXing, bSubEvents, eSubEvents).isSuccess()) &&
323 hitCollList.empty()) {
324 ATH_MSG_ERROR("Could not fill TimedHitCollList");
325 return StatusCode::FAILURE;
326 } else {
327 ATH_MSG_VERBOSE(hitCollList.size() << " RPCSimHitCollection with key " << m_inputHitCollectionName << " found");
328 }
329
330 TimedHitCollList::iterator iColl(hitCollList.begin());
331 TimedHitCollList::iterator endColl(hitCollList.end());
332
333 // Iterating over the list of collections
334 for (; iColl != endColl; ++iColl) {
335 RPCSimHitCollection* hitCollPtr = new RPCSimHitCollection(*iColl->second);
336 PileUpTimeEventIndex timeIndex(iColl->first);
337
338 ATH_MSG_DEBUG("RPCSimHitCollection found with " << hitCollPtr->size() << " hits");
339 ATH_MSG_VERBOSE("time index info. time: " << timeIndex.time() << " index: " << timeIndex.index() << " type: " << timeIndex.type());
340
341 m_thpcRPC->insert(timeIndex, hitCollPtr);
342 m_RPCHitCollList.emplace_back(hitCollPtr);
343 }
344
345 return StatusCode::SUCCESS;
346}
347
348//--------------------------------------------
349// Get next event and extract collection of hit collections:
350StatusCode RpcDigitizationTool::getNextEvent(const EventContext& ctx) {
351 ATH_MSG_DEBUG("RpcDigitizationTool::getNextEvent()");
352
353 // initialize pointer
354 m_thpcRPC.reset();
355
356 // get the container(s)
357 using TimedHitCollList = PileUpMergeSvc::TimedList<RPCSimHitCollection>::type;
358
359 // In case of single hits container just load the collection using read handles
360 if (!m_onlyUseContainerName) {
361 SG::ReadHandle<RPCSimHitCollection> hitCollection(m_hitsContainerKey, ctx);
362 if (!hitCollection.isValid()) {
363 ATH_MSG_ERROR("Could not get RPCSimHitCollection container " << hitCollection.name() << " from store "
364 << hitCollection.store());
365 return StatusCode::FAILURE;
366 }
367
368 // create a new hits collection
369 m_thpcRPC = std::make_unique<TimedHitCollection<RPCSimHit>>(1);
370 m_thpcRPC->insert(0, hitCollection.cptr());
371 ATH_MSG_DEBUG("RPCSimHitCollection found with " << hitCollection->size() << " hits");
372
373 return StatusCode::SUCCESS;
374 }
375 // this is a list<pair<time_t, DataLink<RPCSimHitCollection> > >
376 TimedHitCollList hitCollList;
377
378 if (!(m_mergeSvc->retrieveSubEvtsData(m_inputHitCollectionName, hitCollList).isSuccess())) {
379 ATH_MSG_ERROR("Could not fill TimedHitCollList");
380 return StatusCode::FAILURE;
381 }
382 if (hitCollList.empty()) {
383 ATH_MSG_ERROR("TimedHitCollList has size 0");
384 return StatusCode::FAILURE;
385 } else {
386 ATH_MSG_DEBUG(hitCollList.size() << " RPCSimHitCollections with key " << m_inputHitCollectionName << " found");
387 }
388
389 // create a new hits collection
390 m_thpcRPC = std::make_unique<TimedHitCollection<RPCSimHit>>();
391 // now merge all collections into one
392 TimedHitCollList::iterator iColl(hitCollList.begin());
393 TimedHitCollList::iterator endColl(hitCollList.end());
394 while (iColl != endColl) {
395 const RPCSimHitCollection* p_collection(iColl->second);
396 m_thpcRPC->insert(iColl->first, p_collection);
397 // if ( m_debug ) ATH_MSG_DEBUG ( "RPCSimHitCollection found with "
398 // << p_collection->size() << " hits" ); // loop on the hit collections
399 ++iColl;
400 }
401 return StatusCode::SUCCESS;
402}
403
404//--------------------------------------------
405StatusCode RpcDigitizationTool::mergeEvent(const EventContext& ctx) {
406 StatusCode status = StatusCode::SUCCESS;
407
408 ATH_MSG_DEBUG("RpcDigitizationTool::in mergeEvent()");
409 // create and record the Digit container in StoreGate
410 SG::WriteHandle<RpcDigitContainer> digitContainer(m_outputDigitCollectionKey, ctx);
411 ATH_CHECK(digitContainer.record(std::make_unique<RpcDigitContainer>(m_idHelper->module_hash_max())));
412 ATH_MSG_DEBUG("RpcDigitContainer recorded in StoreGate.");
413
414 // Create and record the SDO container in StoreGate
415 SG::WriteHandle<MuonSimDataCollection> sdoContainer(m_outputSDO_CollectionKey, ctx);
416 ATH_CHECK(sdoContainer.record(std::make_unique<MuonSimDataCollection>()));
417 ATH_MSG_DEBUG("RpcSDOCollection recorded in StoreGate.");
418
420 m_sdo_tmp_map.clear();
422
423 Collections_t collections;
424 status = doDigitization(ctx, collections, sdoContainer.ptr());
425 if (status.isFailure()) { ATH_MSG_ERROR("doDigitization Failed"); }
426 for (size_t coll_hash = 0; coll_hash < collections.size(); ++coll_hash) {
427 if (collections[coll_hash]) {
428 ATH_CHECK( digitContainer->addCollection (collections[coll_hash].release(), coll_hash) );
429 }
430 }
431
432 // Clean-up
433 m_RPCHitCollList.clear();
434
435 return status;
436}
437
438//--------------------------------------------
439StatusCode RpcDigitizationTool::processAllSubEvents(const EventContext& ctx) {
440 StatusCode status = StatusCode::SUCCESS;
441
442 // merging of the hit collection in getNextEvent method
443
444 ATH_MSG_DEBUG("RpcDigitizationTool::in digitize()");
445
446 // create and record the Digit container in StoreGate
447 SG::WriteHandle<RpcDigitContainer> digitContainer(m_outputDigitCollectionKey, ctx);
448 ATH_CHECK(digitContainer.record(std::make_unique<RpcDigitContainer>(m_idHelper->module_hash_max())));
449 ATH_MSG_DEBUG("RpcDigitContainer recorded in StoreGate.");
450
451 // Create and record the SDO container in StoreGate
452 SG::WriteHandle<MuonSimDataCollection> sdoContainer(m_outputSDO_CollectionKey, ctx);
453 ATH_CHECK(sdoContainer.record(std::make_unique<MuonSimDataCollection>()));
454 ATH_MSG_DEBUG("RpcSDOCollection recorded in StoreGate.");
455
457 m_sdo_tmp_map.clear();
459
460 if (!m_thpcRPC) {
461 status = getNextEvent(ctx);
462 if (StatusCode::FAILURE == status) {
463 ATH_MSG_INFO("There are no RPC hits in this event");
464 return status; // there are no hits in this event
465 }
466 }
467
468 Collections_t collections;
469 ATH_CHECK(doDigitization(ctx, collections, sdoContainer.ptr()));
470 for (size_t coll_hash = 0; coll_hash < collections.size(); ++coll_hash) {
471 if (collections[coll_hash]) {
472 ATH_CHECK( digitContainer->addCollection (collections[coll_hash].release(), coll_hash) );
473 }
474 }
475
476 return status;
477}
478
479//--------------------------------------------
480StatusCode RpcDigitizationTool::doDigitization(const EventContext& ctx,
481 Collections_t& collections,
482 MuonSimDataCollection* sdoContainer) {
483 ATHRNG::RNGWrapper* rngWrapper = m_rndmSvc->getEngine(this);
484 rngWrapper->setSeed(name(), ctx);
485 CLHEP::HepRandomEngine* rndmEngine = rngWrapper->getEngine(ctx);
486
487 const MuonGM::MuonDetectorManager* detMgr{nullptr};
488 ATH_CHECK(retrieveCondData(ctx, m_detMgrKey, detMgr));
489
490
491 std::unique_ptr<RPCSimHitCollection> inputSimHitColl{std::make_unique<RPCSimHitCollection>("RPC_Hits")};
492
493
494 // get the iterator pairs for this DetEl
495 // iterate over hits
496 TimedHitCollection<RPCSimHit>::const_iterator i, e;
497
498 // Perform null check on m_thpcRPC
499 if (!m_thpcRPC) {
500 ATH_MSG_ERROR("m_thpcRPC is null");
501 return StatusCode::FAILURE;
502 }
503
504 struct SimDataContent {
505 Identifier channelId{};
506 std::vector<MuonSimData::Deposit> deposits;
507 Amg::Vector3D gpos{Amg::Vector3D::Zero()};
508 double simTime{0.};
509 };
510
511 while (m_thpcRPC->nextDetectorElement(i, e)) {
512 // to store the a single
513
514 std::map<Identifier, SimDataContent> channelSimDataMap;
515
516 // Loop over the hits:
517 while (i != e) {
518 ATH_MSG_DEBUG("RpcDigitizationTool::loop over the hits");
519
520 TimedHitPtr<RPCSimHit> phit(*i++);
521
522 // the hit
523 const RPCSimHit& hit(*phit);
524 // the hit id
525 const int idHit = hit.RPCid();
526 // the global time (G4 time + bunch time)
527 const double globalHitTime{hitTime(phit)};
528 // the G4 time or TOF from IP
529 const double G4Time{hit.globalTime()};
530 // the bunch time
531 const double bunchTime{globalHitTime - hit.globalTime()};
532
533 ATH_MSG_DEBUG("Global time " << globalHitTime << " G4 time " << G4Time << " Bunch time " << bunchTime);
534
535 if (!m_simHitValidKey.empty()) {
536 ATH_MSG_VERBOSE("Validation: globalHitTime, G4Time, BCtime = " << globalHitTime << " " << G4Time << " " << bunchTime);
537 inputSimHitColl->Emplace(idHit, globalHitTime, hit.localPosition(),
538 HepMcParticleLink::getRedirectedLink(phit->particleLink(), phit.eventId(), ctx), // This link should now correctly resolve to the TruthEvent McEventCollection in the main StoreGateSvc.
539 hit.postLocalPosition(),
540 hit.energyDeposit(), hit.stepLength(), hit.particleEncoding(), hit.kineticEnergy());
541 }
542
543 // convert sim id helper to offline id
544 const std::string stationName = m_muonHelper->GetStationName(idHit);
545 const int stationEta = m_muonHelper->GetZSector(idHit);
546 const int stationPhi = m_muonHelper->GetPhiSector(idHit);
547 const int doubletR = m_muonHelper->GetDoubletR(idHit);
548 const int doubletZ = m_muonHelper->GetDoubletZ(idHit);
549 const int doubletPhi = m_muonHelper->GetDoubletPhi(idHit);
550 int gasGap = m_muonHelper->GetGasGapLayer(idHit);
551
552 if (m_muonHelper->GetMeasuresPhi(idHit)) continue; // Skip phi strip . To be created after efficiency evaluation
553
554
555 bool isValid{false};
556 const Identifier elementID = m_idHelper->elementID(stationName,stationEta,stationPhi,doubletR, isValid);
557 if (!isValid) {
558 ATH_MSG_WARNING("Failed to construct the element ID from "<<stationName
559 <<", stationEta: "<<stationEta<<", stationPhi: "<<stationPhi<<", doubletR: "<<doubletR);
560 continue;
561 }
562 // construct Atlas identifier from components
563 ATH_MSG_DEBUG("creating id for hit in element:"
564 << " stationName " << stationName << " stationEta " << stationEta << " stationPhi " << stationPhi << " doubletR "
565 << doubletR << " doubletZ " << doubletZ << " doubletPhi " << doubletPhi << " gasGap " << gasGap);
566 const Identifier detElId{m_idHelper->channelID(elementID, doubletZ, doubletPhi, 1,0, 1, isValid)};
567 if (!isValid) {
568 continue;
569 }
570 const RpcReadoutElement* reEle = detMgr->getRpcReadoutElement(detElId);
572 if (false && reEle->rotatedRpcModule()) {
573 gasGap = gasGap == 1 ? 2 : 1;
574 }
575
576
577 bool isValidEta{false}, isValidPhi{false};
578 const Identifier idpaneleta = m_idHelper->channelID(elementID, doubletZ, doubletPhi, gasGap, 0, 1, isValidEta);
579 const Identifier idpanelphi = m_idHelper->channelID(elementID, doubletZ, doubletPhi, gasGap, 1, 1, isValidPhi);
580 if (!isValidEta || !isValidPhi) {
581 ATH_MSG_WARNING("Found an invalid identifier "
582 << " stationName " << stationName << " stationEta " << stationEta << " stationPhi " << stationPhi
583 << " doubletR " << doubletR << " doubletZ " << doubletZ << " doubletPhi " << doubletPhi << " gasGap "
584 << gasGap);
585 continue;
586 }
587 // loop on eta and phi to apply correlated efficiency between the two views
588
590 const double tmp_CorrJitter = m_idHelper->stationName(idpaneleta) < 2 ? m_CorrJitter_BIS78 : m_CorrJitter;
592 const double corrtimejitter = tmp_CorrJitter > 0.01 ?
593 CLHEP::RandGaussZiggurat::shoot(rndmEngine, 0., tmp_CorrJitter) : 0.; // correlated jitter
594 // handle here the special case where eta panel is dead => phi strip status (dead or eff.) cannot be resolved;
595 // measured panel eff. will be used in that case and no phi strip killing will happen
596
597
598 // Extrapolate the hit to the gas gap centre located at x=0
599 const Amg::Vector3D hitDir{(hit.postLocalPosition() - hit.localPosition()).unit()};
600 const Amg::Vector3D gapCentre = hit.localPosition() +
601 Amg::intersect<3>(hit.localPosition(), hitDir, Amg::Vector3D::UnitX(), 0).value_or(0) * hitDir;
602
603 std::array<int, 3> pcseta = physicalClusterSize(ctx, reEle, idpaneleta, gapCentre, rndmEngine); // set to one for new algorithms
604 ATH_MSG_VERBOSE("Simulated cluster on eta panel: size/first/last= " << pcseta[0] << "/" << pcseta[1] << "/" << pcseta[2]);
605 std::array<int, 3> pcsphi = physicalClusterSize(ctx, reEle, idpanelphi, gapCentre, rndmEngine); // set to one for new algorithms
606 ATH_MSG_VERBOSE("Simulated cluster on phi panel: size/first/last= " << pcsphi[0] << "/" << pcsphi[1] << "/" << pcsphi[2]);
607
608
609
610 // create Identifiers
611 const Identifier atlasRpcIdeta = m_idHelper->channelID(elementID, doubletZ, doubletPhi, gasGap, 0, pcseta[1], isValidEta);
612 const Identifier atlasRpcIdphi = m_idHelper->channelID(elementID, doubletZ, doubletPhi, gasGap, 1, pcsphi[1], isValidPhi);
613
614 const HepMcParticleLink particleLink = HepMcParticleLink::getRedirectedLink(phit->particleLink(), phit.eventId(), ctx); // This link should now correctly resolve to the TruthEvent McEventCollection in the main StoreGateSvc.
615 const auto [etaStripOn, phiStripOn] = detectionEfficiency(ctx, idpaneleta, idpanelphi, rndmEngine, particleLink);
616 ATH_MSG_DEBUG("SetPhiOn " << phiStripOn << " SetEtaOn " << etaStripOn);
617
618 for (bool imeasphi : {false, true}) {
619 if (!imeasphi && (!etaStripOn || !isValidEta)) continue;
620 if (imeasphi && (!phiStripOn || !isValidPhi)) continue;
621
622
623 // get Identifier and list of clusters for this projection
624 const Identifier& atlasId = !imeasphi ? atlasRpcIdeta : atlasRpcIdphi;
625 std::array<int, 3> pcs{!imeasphi ? pcseta : pcsphi};
626
627 ATH_MSG_DEBUG("SetOn: stationName " << stationName << " stationEta " << stationEta << " stationPhi " << stationPhi
628 << " doubletR " << doubletR << " doubletZ " << doubletZ << " doubletPhi " << doubletPhi
629 << " gasGap " << gasGap << " measphi " << imeasphi);
630
631 // pcs contains the cluster size, the first strip number and the last strip number of the cluster
632 pcs = TurnOnStrips(reEle, std::move(pcs), atlasId);
633 if (pcs[2] < 0){
634 continue;
635 }
636
637 ATH_MSG_DEBUG("Simulated cluster1: size/first/last= " << pcs[0] << "/" << pcs[1] << "/" << pcs[2]);
638
639
640 const Amg::Vector3D pos = fromSimHitToLayer(reEle, atlasId) * hit.localPosition();
641 const Amg::Vector3D gpos = reEle->transform(atlasId) * pos;
642
643 ATH_MSG_VERBOSE(" evt: "<<ctx.eventID().event_number()
644 <<" hit "<<m_idHelper->print_to_string(atlasId)
645 <<" local simHit "<<Amg::toString(hit.localPosition())
646 <<" corrected: "<<Amg::toString(pos)
647 <<" transform: "<<GeoTrf::toString(fromSimHitToLayer(reEle, atlasId))
648 <<" local strip: "<<Amg::toString(reEle->localToGlobalTransf(atlasId).inverse()*reEle->stripPos(atlasId))
649 <<" local strip (II): "<<Amg::toString(reEle->transform(atlasId).inverse()*reEle->stripPos(atlasId))
650 <<" global: "<<Amg::toString(gpos)
651 <<" strip Pos: "<<Amg::toString(reEle->stripPos(atlasId)));
652
653 // Calculate propagation time along readout strip in seconds
654 double proptime = PropagationTime(reEle, atlasId, gpos);
655
656 double tns = G4Time + proptime + corrtimejitter; // the time is in nanoseconds
657 ATH_MSG_VERBOSE("TOF+propagation time " << tns << " /s where proptime " << proptime << "/s");
658
659 double time = tns + bunchTime;
660 ATH_MSG_VERBOSE("final time in ns: BC+TOF+prop " << time << " /ns");
661
662 // pack propagation time along strip, bunch time and local hit position
663 long long int packedMCword = PackMCTruth(proptime, bunchTime, pos.y(), pos.z());
664 //cppcheck-suppress invalidPointerCast
665 double* b = reinterpret_cast<double*>(&packedMCword);
666
668 // create here deposit for MuonSimData
669 // MuonMCData first word is the packing of : proptime, bunchTime, posy, posz
670 // MuonMCData second word is the total hit time: bunchcTime+tof+proptime+correlatedJitter / ns
671 MuonSimData::Deposit deposit(particleLink, MuonMCData((*b), time)); // store tof+strip_propagation+corr.jitter
672 // MuonMCData((*b),G4Time+bunchTime+proptime )); // store tof+strip_propagation
673
674 // Do not store pile-up truth information
675 if (m_includePileUpTruth || !HepMC::ignoreTruthLink(phit->particleLink(), m_vetoPileUpTruthLinks)) {
676 if (std::abs(hit.particleEncoding()) == 13 || hit.particleEncoding() == 0) {
677 if (channelSimDataMap.find(atlasId) == channelSimDataMap.end()) {
678 SimDataContent& content = channelSimDataMap[atlasId];
679 content.channelId = atlasId;
680 content.deposits.push_back(deposit);
681 content.gpos = reEle->transform(atlasId)*
682 fromSimHitToLayer(reEle,atlasId) * gapCentre;
683 content.simTime = hitTime(phit);
684 ATH_MSG_VERBOSE("adding SDO entry: r " << content.gpos.perp() << " z " << content.gpos.z());
685 }
686 }
687 }
688
689
690 //---------------------------------------------------------------------
691 // construct new digit and store it in the respective digit collection
692 // --------------------------------------------------------------------
693
694 // we create one digit-vector/deposit for each strip in the cluster
695 bool isValid{false};
696 for (int clus = pcs[1]; clus <= pcs[2]; ++clus) {
697 Identifier newId = m_idHelper->channelID(stationName, stationEta, stationPhi, doubletR, doubletZ,
698 doubletPhi, gasGap, imeasphi, clus, isValid);
699 if (!isValid) {
700 ATH_MSG_WARNING(__FILE__<<":"<<__LINE__<< "Channel "<< stationName<<" "<<stationEta<<" "<<stationPhi<<" "<< doubletR<<" "<<doubletZ
701 <<" "<< doubletPhi<<" "<< gasGap <<" "<< imeasphi<<" "<< clus<<" is invalid");
702 continue;
703 }
704
705 if (!m_idHelper->valid(newId)) {
706 if (stationName.find("BI") != std::string::npos) {
707 ATH_MSG_WARNING("Temporary skipping creation of RPC digit for stationName="
708 << stationName << ", eta=" << stationEta << ", phi=" << stationPhi << ", doubletR=" << doubletR
709 << ", doubletZ=" << doubletZ << ", doubletPhi=" << doubletPhi << ", gasGap=" << gasGap
710 << ", measuresPhi=" << imeasphi << ", strip=" << clus << ", cf. ATLASRECTS-6124");
711 return StatusCode::SUCCESS;
712 } else {
713 ATH_MSG_ERROR("Created an invalid id, aborting!");
714 m_idHelper->print(newId);
715 return StatusCode::FAILURE;
716 }
717 }
718
722 // One identifier but several deposits // name m_sdo_tmp_map is wrong call it m_sdo_map
723 if (m_sdo_tmp_map.find(newId) == m_sdo_tmp_map.end()) {
724 std::vector<MuonSimData::Deposit> newdeps;
725 newdeps.push_back(deposit);
726 m_sdo_tmp_map.insert(std::map<Identifier, std::vector<MuonSimData::Deposit>>::value_type(newId, newdeps));
727 } else {
728 m_sdo_tmp_map[newId].push_back(deposit);
729 }
730 } // end for cluster
731 } // loop on eta and phi
732 } // end loop hits
733
734 if (m_muonOnlySDOs) {
735 for (auto it = channelSimDataMap.begin(); it != channelSimDataMap.end(); ++it) {
736 MuonSimData simData(it->second.deposits, 0);
737 simData.setPosition(it->second.gpos);
738 simData.setTime(it->second.simTime);
739 auto insertResult = sdoContainer->insert(std::make_pair(it->first, simData));
740 if (!insertResult.second)
741 ATH_MSG_WARNING("Attention: this sdo is not recorded, since the identifier already exists in the sdoContainer map");
742 }
743 }
744
745 } // end loop detector elements
746
748
749 std::map<Identifier, std::vector<MuonSimData::Deposit>>::iterator map_iter = m_sdo_tmp_map.begin();
750 ATH_MSG_DEBUG("Start the digit map loop");
751
752 for (; map_iter != m_sdo_tmp_map.end(); ++map_iter) {
753 // Identifier
754 const Identifier theId = (*map_iter).first;
755 ATH_MSG_DEBUG("in the map loop: id " << m_idHelper->show_to_string(theId));
756 // Deposit
757 const std::vector<MuonSimData::Deposit> theDeps = (*map_iter).second;
758
759 // store the SDO from the muon
760 MuonSimData::Deposit theMuon; // useful beacuse it sorts the digits in ascending time.
761 std::multimap<double, MuonSimData::Deposit> times; // extract here time info from deposits.
762
763 // loop on the vector deposit
764 for (unsigned int k = 0; k < theDeps.size(); k++) {
765 double time = theDeps[k].second.secondEntry();
766 times.insert(std::multimap<double, MuonSimData::Deposit>::value_type(time, theDeps[k]));
767 }
768
769 // now iterate again over the multimap entries and store digits after dead time applied
770
771 IdContext rpcContext = m_idHelper->module_context(); // work on chamber context
772
773 std::multimap<double, MuonSimData::Deposit>::iterator map_dep_iter = times.begin();
774
775 // loop to suppress digits too close in time (emulate Front-End and CMA dead time)
776 double last_time = -10000; // init to high value
777 for (; map_dep_iter != times.end(); ++map_dep_iter) {
778 double currTime = (*map_dep_iter).first;
779 ATH_MSG_VERBOSE("deposit with time " << currTime);
780
781 if (!m_muonOnlySDOs && !m_sdoAreOnlyDigits) {
782 // store (before any cut: all G4 hits) in the SDO container
783 // Identifier sdo and digit are the same
784 if (sdoContainer->find(theId) != sdoContainer->end()) // Identifier exist -> increase deposit
785 {
786 std::map<Identifier, MuonSimData>::const_iterator it = sdoContainer->find(theId);
787 std::vector<MuonSimData::Deposit> deps = ((*it).second).getdeposits();
788 deps.push_back((*map_dep_iter).second);
789 } else // Identifier does not exist -> create (Id,deposit)
790 {
791 std::vector<MuonSimData::Deposit> deposits;
792 deposits.push_back((*map_dep_iter).second);
793 std::pair<std::map<Identifier, MuonSimData>::iterator, bool> insertResult =
794 sdoContainer->insert(std::make_pair(theId, MuonSimData(deposits, 0)));
795 if (!insertResult.second)
796 ATH_MSG_ERROR(
797 "Attention TEMP: this sdo is not recorded, since the identifier already exists in the sdoContainer map");
798 }
799 }
800 // apply dead time
801 if (std::abs(currTime - last_time) > (m_deadTime)) {
802 ATH_MSG_DEBUG("deposit with time " << currTime << " is distant enough from previous (if any) hit on teh same strip");
803 last_time = (*map_dep_iter).first;
804
805 // first add time jitter to the time:
806 double uncorrjitter = 0;
807 double tmp_UncorrJitter = m_UncorrJitter;
808 if (m_idHelper->stationName(theId) < 2) tmp_UncorrJitter = m_UncorrJitter_BIS78;
809 if (tmp_UncorrJitter > 0.01) uncorrjitter = CLHEP::RandGaussZiggurat::shoot(rndmEngine, 0., tmp_UncorrJitter);
810 // Historically patch for the cavern background
811 // Now we subtract TOF from IP to assume full time calibrated detector (t=0 for particle from IP at light speed)
812 // We add a time shift to emulate FE global offset
813
814 const RpcReadoutElement* ele = detMgr->getRpcReadoutElement(theId);
815 Amg::Vector3D posi = ele->stripPos(theId);
816 double tp = m_patch_for_rpc_time ? posi.mag() / Gaudi::Units::c_light : 0.;
817 // Calculate propagation time for a hit at the center of the strip, to be subtructed as well as the nominal TOF
818 double propTimeFromStripCenter = PropagationTime(ele, theId, posi);
819 double newDigit_time = currTime + uncorrjitter + m_rpc_time_shift - tp - propTimeFromStripCenter;
820
821 double digi_ToT = -1.; // Time over threshold, for Narrow-gap RPCs only
822 if (m_idHelper->stationName(theId) < 2) digi_ToT = timeOverThreshold(rndmEngine); //mn
823
824 ATH_MSG_VERBOSE("last_time=currTime " << last_time << " jitter " << uncorrjitter << " TOFcorrection " << tp << " shift "
825 << m_rpc_time_shift << " newDigit_time " << newDigit_time);
826
827 // Apply readout window (sensitive detector time window)
828 bool outsideDigitizationWindow = outsideWindow(newDigit_time);
829 if (outsideDigitizationWindow) {
830 ATH_MSG_VERBOSE("hit outside digitization window - do not produce digits");
831 ATH_MSG_DEBUG("Hit outside time window!!"
832 << " hit time (ns) = " << newDigit_time << " timeWindow = " << m_timeWindowLowerOffset << " / "
833 << m_timeWindowUpperOffset);
834
835 continue;
836 }
837 // ok, let's store this digit
838 // this is an accepted hit to become digit
839 last_time = (*map_dep_iter).first;
840
841 std::unique_ptr<RpcDigit> newDigit = std::make_unique<RpcDigit>(theId, newDigit_time, digi_ToT, false);
842
843 Identifier elemId = m_idHelper->elementID(theId);
844 RpcDigitCollection* digitCollection = nullptr;
845
846 IdentifierHash coll_hash;
847 if (m_idHelper->get_hash(elemId, coll_hash, &rpcContext)) {
848 ATH_MSG_ERROR("Unable to get RPC hash id from RPC Digit collection "
849 << "context begin_index = " << rpcContext.begin_index()
850 << " context end_index = " << rpcContext.end_index() << " the identifier is ");
851 elemId.show();
852 }
853
854 // make new digit
855 ATH_MSG_DEBUG("Digit Id = " << m_idHelper->show_to_string(theId) << " digit time " << newDigit_time);
856
857 // remember new collection.
858 if (coll_hash >= collections.size()) {
859 collections.resize (coll_hash+1);
860 }
861 digitCollection = collections[coll_hash].get();
862 if (!digitCollection) {
863 collections[coll_hash] = std::make_unique<RpcDigitCollection>(elemId, coll_hash);
864 digitCollection = collections[coll_hash].get();
865 }
866 digitCollection->push_back(std::move(newDigit));
867
868 if (!m_muonOnlySDOs && m_sdoAreOnlyDigits) {
869 // put SDO collection in StoreGate
870 if (sdoContainer->find(theId) != sdoContainer->end()) {
871 std::map<Identifier, MuonSimData>::const_iterator it = sdoContainer->find(theId);
872 std::vector<MuonSimData::Deposit> deps = ((*it).second).getdeposits();
873 deps.push_back((*map_dep_iter).second);
874 } else {
875 std::vector<MuonSimData::Deposit> deposits;
876 deposits.push_back((*map_dep_iter).second);
877 std::pair<std::map<Identifier, MuonSimData>::iterator, bool> insertResult =
878 sdoContainer->insert(std::make_pair(theId, MuonSimData(deposits, 0)));
879 if (!insertResult.second)
880 ATH_MSG_ERROR(
881 "Attention: this sdo is not recorded, since teh identifier already exists in the sdoContainer map");
882 }
883 }
884
885 } else
886 ATH_MSG_DEBUG("discarding digit due to dead time: " << (*map_dep_iter).first << " " << last_time);
887 }
888
889 } // loop to suppress digits too close in time ended
890
891 // reset the pointer if it not null
892 m_thpcRPC.reset();
893 if (!m_simHitValidKey.empty()) {
894 SG::WriteHandle<RPCSimHitCollection> validHandle{m_simHitValidKey, ctx};
895 ATH_CHECK(validHandle.record(std::move(inputSimHitColl)));
896 }
897
898 return StatusCode::SUCCESS;
899}
900Amg::Transform3D RpcDigitizationTool::fromSimHitToLayer(const MuonGM::RpcReadoutElement* reEle,
901 const Identifier& layerId) const {
902
903 Amg::Vector3D lGasGapPos = reEle->localGasGapPos(layerId);
904 if (reEle->NphiStripPanels() != reEle->nGasGapPerLay()) {
905 lGasGapPos.y() =0.;
906 }
907
911 const bool flip = reEle->numberOfLayers() == 2 &&
912 (m_idHelper->gasGap(layerId) == 2) != reEle->rotatedRpcModule();
913 const Amg::Transform3D fromHitToGap{reEle->transform(layerId).inverse() *
914 reEle->absTransform() * Amg::getTranslate3D(lGasGapPos) *
915 (flip ? Amg::getRotateY3D(180.*Gaudi::Units::deg) : Amg::Transform3D::Identity())};
916 ATH_MSG_VERBOSE("Transformation to go from hit to gap restframe "<<m_idHelper->print_to_string(layerId)
917 <<" "<<Amg::toString(fromHitToGap));
918 return fromHitToGap;
919}
920
921//--------------------------------------------
922std::array<int, 3> RpcDigitizationTool::physicalClusterSize(const EventContext& ctx,
923 const RpcReadoutElement* ele,
924 const Identifier& id,
925 const Amg::Vector3D& gapCentre,
926 CLHEP::HepRandomEngine* rndmEngine) const {
927
928
929 std::array<int, 3> result{};
930
931 const Amg::Vector3D position = fromSimHitToLayer(ele, id) * gapCentre;
932
933 const int doubletPhi = m_idHelper->doubletPhi(id);
934 const int gasGap = m_idHelper->gasGap(id);
935 const bool measuresPhi = m_idHelper->measuresPhi(id);
936 const double pitch= ele->StripPitch(measuresPhi);
937
938
939 const int nstrip = ele->stripNumber(position.block<2,1>(0,0), id);
940 const int numStrips = ele->Nstrips(measuresPhi);
941
942 result[1] = nstrip;
943 result[2] = nstrip;
944
945 if (nstrip < 1 || nstrip > numStrips) {
946 return make_array<int, 3>(-1);
947 }
948 const Amg::Vector3D locStripPos = ele->transform(id).inverse()*ele->stripPos(doubletPhi, gasGap, measuresPhi, nstrip);
949 float xstripnorm = (locStripPos -position).x() / pitch ;
950 result[0] = determineClusterSize(ctx, id, xstripnorm, rndmEngine);
951
952 //
953
954
955 if (m_turnON_clustersize == false) result[0] = 1;
956
957 return result;
958}
959
960//--------------------------------------------
961std::array<int, 3> RpcDigitizationTool::TurnOnStrips(const RpcReadoutElement* ele,
962 std::array<int, 3>&& pcs,
963 const Identifier& id) const {
964
965
966 const int nstrips = ele->Nstrips(m_idHelper->measuresPhi(id));
967
968 if (pcs[0] == -2) {
969 pcs[1] = pcs[2] - 1;
970 } else if (pcs[0] == 2) {
971 pcs[2] = pcs[1] + 1;
972 } else if (pcs[0] > 2) {
973 pcs[1] = pcs[1] - pcs[0] / 2;
974 if (fmod(pcs[0], 2) == 0) pcs[1] = pcs[1] + 1;
975 pcs[2] = pcs[1] + pcs[0] - 1;
976 } else if (pcs[0] < -2) {
977 pcs[1] = pcs[1] + pcs[0] / 2;
978 pcs[2] = pcs[1] - pcs[0] - 1;
979 }
980
981 // cut the clusters at the beginning and at the end of the chamber
982
983 pcs[1] = std::clamp(pcs[1], 1, nstrips);
984 pcs[2] = std::clamp(pcs[2], 1, nstrips);
985
986 pcs[0] = pcs[2] - pcs[1] + 1;
987
988 return pcs;
989}
990
991//--------------------------------------------
992double RpcDigitizationTool::PropagationTime(const MuonGM::RpcReadoutElement* ele,
993 const Identifier& id,
994 const Amg::Vector3D& globPos) const {
995
996 double distance{0.};
997 if (m_idHelper->measuresPhi(id)) {
998 distance = ele->distanceToPhiReadout(globPos);
999 } else {
1000 distance = ele->distanceToEtaReadout(globPos);
1001 }
1002
1003 // distance in mm, SIG_VEL in ns/m
1004 return std::abs(distance * SIG_VEL * 1.e-3);
1005}
1006
1007//--------------------------------------------
1008long long int RpcDigitizationTool::PackMCTruth(float proptime, float bctime, float posy, float posz) const {
1009 // start with proptime: it is usually ~ns. It comes in ns. We express it in ns/10. use only 8 bits
1010 if (proptime < 0) {
1011 ATH_MSG_WARNING("A poblem: packing a propagation time <0 " << proptime << " redefine it as 0");
1012 proptime = 0.;
1013 }
1014 long long int new_proptime = int(proptime * 10) & 0xff;
1015
1016 // now tof. it is ~100ns. comes in ns. express it in ns/10. 16 bits needed (0-32768)
1017 // now BC time: it is ~100ns. comes in ns. express it in ns/10. 16 bits needed (0-32768)
1018 // can be negative (=> add 300 ns)
1019
1020 long long int new_bctime = int((bctime + 300.) * 10.) & 0xffff;
1021
1022 // posy: ~1000mm comes in mm, write it in mm*10. need 16 bits (0-32768)
1023 // can be negative (=>add 1500 mm)
1024
1025 long long int new_posy = int((posy + 1500.) * 10.) & 0xffff;
1026
1027 // posz: ~1000mm comes in mm, write it in mm*10. need 16 bits (0-32768)
1028 // can be negative (=>add 1500 mm)
1029
1030 long long int new_posz = int((posz + 1500.) * 10.) & 0xffff;
1031
1032 return (new_proptime + (new_bctime << 8) + (new_posy << 24) + (new_posz << 40));
1033}
1034
1035//--------------------------------------------
1036void RpcDigitizationTool::UnPackMCTruth(double theWord, float& proptime, float& bctime, float& posy, float& posz) {
1037 // int64_t is just a shorter way of writing long long int
1038 using Repacker = union
1039
1040 {
1041 double dWord;
1042
1043 int64_t iWord;
1044 };
1045 Repacker MCTruth;
1046 MCTruth.dWord = theWord;
1047 proptime = ((MCTruth.iWord) & 0x00000000000000ffLL) / 10.;
1048 bctime = (((MCTruth.iWord) & 0x0000000000ffff00LL) >> 8) / 10.;
1049 posy = (((MCTruth.iWord) & 0x000000ffff000000LL) >> 24) / 10.;
1050 posz = (((MCTruth.iWord) & 0x00ffff0000000000LL) >> 40) / 10.;
1051
1052 //
1053 bctime = bctime - 300.;
1054 posy = posy - 1500.;
1055 posz = posz - 1500.;
1056}
1057
1058//--------------------------------------------
1059StatusCode RpcDigitizationTool::fillTagInfo() {
1060 // get TagInfoMgr
1061 SmartIF<ITagInfoMgr> tagInfoMgr{Gaudi::svcLocator()->service("TagInfoMgr")}; // Tag Info Manager
1062 if (!tagInfoMgr) { return StatusCode::FAILURE; }
1063
1064 std::string RpctimeSchema = "";
1065 std::stringstream RpctimeShift;
1066 RpctimeShift << (int)m_rpc_time_shift;
1067
1068 if (m_patch_for_rpc_time) {
1069 RpctimeSchema = "Datalike_TOFoff_TimeShift" + RpctimeShift.str() + "nsec";
1070 } else {
1071 RpctimeSchema = "G4like_TOFon_TimeShift" + RpctimeShift.str() + "nsec";
1072 }
1073
1074 StatusCode sc = tagInfoMgr->addTag(m_RPC_TimeSchema, RpctimeSchema);
1075
1076 if (sc.isFailure()) {
1077 ATH_MSG_WARNING(m_RPC_TimeSchema << " " << RpctimeSchema << " not added to TagInfo ");
1078 return sc;
1079 }
1080
1081 ATH_MSG_DEBUG(m_RPC_TimeSchema << " " << RpctimeSchema << " added to TagInfo ");
1082
1083 return StatusCode::SUCCESS;
1084}
1085
1086
1087//--------------------------------------------
1088std::pair<bool,bool> RpcDigitizationTool::detectionEfficiency(const EventContext& ctx,
1089 const Identifier& IdEta,
1090 const Identifier& IdPhi,
1091 CLHEP::HepRandomEngine* rndmEngine,
1092 const HepMcParticleLink& trkParticle) const {
1093
1094
1095
1096 ATH_MSG_DEBUG("RpcDigitizationTool::in DetectionEfficiency");
1097
1098 ATH_MSG_DEBUG("EtaPanelId to look for Eff is " << m_idHelper->show_to_string(IdEta));
1099 ATH_MSG_DEBUG("PhiPanelId to look for Eff is " << m_idHelper->show_to_string(IdPhi));
1100
1101
1102 // dead spacers are not simulated in GEANT4 => their effect must be emulated in the digitizer as an effective max. efficiency = 99%
1103 // (spacers are 1x1cm^2 over a grid of 10x10cm^2 =? geometrical ineff. introduced is 1% for normal incidence)
1104 float maxGeomEff{0.99}, PhiAndEtaEff{0.99}, OnlyEtaEff{0.f}, OnlyPhiEff{0.f};
1105
1106 // 2=BML,3=BMS,4=BOL,5=BOS,8=BMF,9=BOF,10=BOG
1107 int stationName = m_idHelper->stationName(IdEta);
1108 int stationEta = m_idHelper->stationEta(IdEta);
1109 int doubletR = m_idHelper->doubletR(IdEta);
1110
1111 // remove feet extension. driven by joboption
1112 if (m_BOG_BOF_DoubletR2_OFF && (stationName == m_BOF_id || stationName == m_BOG_id) && doubletR == 2) {
1113 return std::make_pair(false, false);
1114 }
1115
1116
1117 if (!m_turnON_efficiency) {
1118 return std::make_pair(true, true);
1119 }
1120 bool etaStripOn{true}, phiStripOn{true};
1121
1122 // int stripetadead = 0 ; // not used
1123 // int stripphidead = 0 ; // not used
1124
1125 unsigned int index = stationName - 2;
1126 // BML and BMS, BOL and BOS come first (stationName= 2 and 3, 4 and 5 -> index 0-3)
1127 if (stationName > 5 && stationName < 50) index = index - 2;
1128 // BMF, BOF and BOG are 8,9,10 => must be 4,5 and 6
1129 else if (stationName > 50)
1130 index = index - 44;
1131 // BME and BOE 53 and 54 are at indices 7 and 8
1132
1133 if (!m_Efficiency_fromCOOL && stationName >= 2) {
1134 if (index > m_PhiAndEtaEff_A.size() || index > m_OnlyEtaEff_A.size() || index > m_OnlyPhiEff_A.size()) {
1135 THROW_EXCEPTION("Index out of array in Detection Efficiency SideA " << index << " stationName = " << stationName);
1136 }
1137
1138 PhiAndEtaEff = m_PhiAndEtaEff_A[index];
1139 OnlyEtaEff = m_OnlyEtaEff_A[index];
1140 OnlyPhiEff = m_OnlyPhiEff_A[index];
1141
1142 if (stationEta < 0) {
1143 if (index > m_PhiAndEtaEff_C.size() || index > m_OnlyEtaEff_C.size() || index > m_OnlyPhiEff_C.size()) {
1144 THROW_EXCEPTION("Index out of array in Detection Efficiency SideC " << index << " stationName = " << stationName);
1145 }
1146 PhiAndEtaEff = m_PhiAndEtaEff_C[index];
1147 OnlyEtaEff = m_OnlyEtaEff_C[index];
1148 OnlyPhiEff = m_OnlyPhiEff_C[index];
1149 }
1150 } else if (stationName < 2 && (!m_Efficiency_fromCOOL || !m_Efficiency_BIS78_fromCOOL)) { // BIS
1151 PhiAndEtaEff = m_PhiAndEtaEff_BIS78;
1152 OnlyEtaEff = m_OnlyEtaEff_BIS78;
1153 OnlyPhiEff = m_OnlyPhiEff_BIS78;
1154 } else { // Efficiency from Cool
1155
1156 const RpcCondDbData* readCdo{nullptr};
1157 if(!retrieveCondData(ctx, m_readKey, readCdo).isSuccess()){
1158 THROW_EXCEPTION("Failed to retrieve conditions object");
1159 }
1160
1161 ATH_MSG_DEBUG("Efficiencies and cluster size + dead strips will be extracted from COOL");
1162
1163 double FracDeadStripEta{0.}, FracDeadStripPhi{0.};
1164 double EtaPanelEfficiency{1.}, PhiPanelEfficiency{1.}, GapEfficiency{1.};
1165 int RPC_ProjectedTracksEta = 0;
1166
1167 std::optional<double> fracDeadStripEtaFromCOOL = readCdo->getFracDeadStrip(IdEta);
1168 std::optional<double> fracDeadStripPhiFromCOOL = readCdo->getFracDeadStrip(IdPhi);
1169
1170 bool noEntryInDb = !fracDeadStripEtaFromCOOL || !fracDeadStripPhiFromCOOL;
1171
1172 FracDeadStripEta = fracDeadStripEtaFromCOOL.value_or(0.);
1173 FracDeadStripPhi = fracDeadStripPhiFromCOOL.value_or(0.);
1174 RPC_ProjectedTracksEta = readCdo->getProjectedTrack(IdEta).value_or(0);
1175
1176 EtaPanelEfficiency = readCdo->getEfficiency(IdEta).value_or(1.);
1177 PhiPanelEfficiency = readCdo->getEfficiency(IdPhi).value_or(1.);
1178 GapEfficiency = readCdo->getGapEfficiency(IdEta).value_or(1.);
1179
1180 if (std::abs(FracDeadStripEta - 1.) < 0.001) {
1181 ATH_MSG_DEBUG("Watch out: SPECIAL CASE: Read from Cool: FracDeadStripEta/Phi "
1182 << FracDeadStripEta << "/" << FracDeadStripPhi << " RPC_ProjectedTracksEta " << RPC_ProjectedTracksEta
1183 << " Eta/PhiPanelEfficiency " << EtaPanelEfficiency << "/" << PhiPanelEfficiency << " gapEff " << GapEfficiency
1184 << " for gas gap " << m_idHelper->show_to_string(IdEta) << " id " << IdEta.get_identifier32().get_compact());
1185 // dead eta panel => cannot determine the strip status for phi strips
1186 // FracDeadStripPhi must be reset to 0. and undefinedPhiStripStatus = true
1187 FracDeadStripPhi = 0.;
1188 ATH_MSG_VERBOSE("Watch out: SPECIAL CASE: Resetting FracDeadStripPhi " << FracDeadStripPhi << " ignoring phi dead strips ");
1189 }
1190
1191 // special test
1192 // here redefining the efficiencies:
1193 // EtaPanelEfficiency = 0.92;
1194 // PhiPanelEfficiency = 0.85;
1195 // GapEfficiency = 0.97;
1196 bool changing = false;
1197 ATH_MSG_DEBUG("Read from Cool: FracDeadStripEta/Phi " << FracDeadStripEta << "/" << FracDeadStripPhi << " RPC_ProjectedTracksEta "
1198 << RPC_ProjectedTracksEta << " Eta/PhiPanelEfficiency " << EtaPanelEfficiency
1199 << "/" << PhiPanelEfficiency << " gapEff " << GapEfficiency);
1200 // if ((1.-FracDeadStripEta)<EtaPanelEfficiency)
1201 if ((maxGeomEff - FracDeadStripEta) - EtaPanelEfficiency < -0.011) {
1202 ATH_MSG_DEBUG("Ineff. from dead strips on Eta Panel larger that measured efficiency: deadFrac="
1203 << FracDeadStripEta << " Panel Eff=" << EtaPanelEfficiency << " for Panel " << m_idHelper->show_to_string(IdEta));
1204 ATH_MSG_DEBUG("... see the corresponding report from RpcDetectorStatusDbTool");
1205 // EtaPanelEfficiency = 1.-FracDeadStripEta;
1206 EtaPanelEfficiency = maxGeomEff - FracDeadStripEta;
1207 changing = true;
1208 }
1209 // if ((1.-FracDeadStripPhi)<PhiPanelEfficiency)
1210 if ((maxGeomEff - FracDeadStripPhi) - PhiPanelEfficiency < -0.011) {
1211 ATH_MSG_DEBUG("Ineff. from dead strips on Phi Panel larger that measured efficiency: deadFrac="
1212 << FracDeadStripPhi << " Panel Eff=" << PhiPanelEfficiency << " for Panel " << m_idHelper->show_to_string(IdPhi));
1213 ATH_MSG_DEBUG("... see the corresponding report among the warnings of RpcDetectorStatusDbTool");
1214 // PhiPanelEfficiency = 1.-FracDeadStripPhi;
1215 PhiPanelEfficiency = maxGeomEff - FracDeadStripPhi;
1216 changing = true;
1217 }
1218 // if ((1.-FracDeadStripEta*FracDeadStripPhi)<GapEfficiency)
1219 if ((maxGeomEff - FracDeadStripEta * FracDeadStripPhi) - GapEfficiency < -0.011) {
1220 ATH_MSG_DEBUG("Ineff. from dead strips on Eta/Phi Panels larger that measured EtaORPhi efficiency: deadFrac="
1221 << FracDeadStripEta * FracDeadStripPhi << " EtaORPhi Eff=" << GapEfficiency << " for GasGap "
1222 << m_idHelper->show_to_string(IdEta));
1223 ATH_MSG_DEBUG("... see the corresponding report among the warnings of RpcDetectorStatusDbTool");
1224 // GapEfficiency = 1.-FracDeadStripEta*FracDeadStripPhi;
1225 GapEfficiency = maxGeomEff - FracDeadStripEta * FracDeadStripPhi;
1226 changing = true;
1227 }
1228 if (changing)
1229 ATH_MSG_DEBUG("Rinormalized Values from Cool: FracDeadStripEta/Phi "
1230 << FracDeadStripEta << "/" << FracDeadStripPhi << " RPC_ProjectedTracksEta " << RPC_ProjectedTracksEta
1231 << " Eta/PhiPanelEfficiency " << EtaPanelEfficiency << "/" << PhiPanelEfficiency << " gapEff " << GapEfficiency);
1232
1233 // gabriele //..stefania - if there are dead strips renormalize the eff. to the active area
1234 if (m_kill_deadstrips) {
1235 if ((FracDeadStripEta > 0.0 && FracDeadStripEta < 1.0) || (FracDeadStripPhi > 0.0 && FracDeadStripPhi < 1.0) || (noEntryInDb)) {
1236 EtaPanelEfficiency = EtaPanelEfficiency / (maxGeomEff - FracDeadStripEta);
1237 PhiPanelEfficiency = PhiPanelEfficiency / (maxGeomEff - FracDeadStripPhi);
1238 GapEfficiency = GapEfficiency / (maxGeomEff - FracDeadStripEta * FracDeadStripPhi);
1239
1240 if (EtaPanelEfficiency > maxGeomEff) EtaPanelEfficiency = maxGeomEff;
1241 if (PhiPanelEfficiency > maxGeomEff) PhiPanelEfficiency = maxGeomEff;
1242 if (GapEfficiency > maxGeomEff) GapEfficiency = maxGeomEff;
1243
1244 if (EtaPanelEfficiency > GapEfficiency) GapEfficiency = EtaPanelEfficiency;
1245 if (PhiPanelEfficiency > GapEfficiency) GapEfficiency = PhiPanelEfficiency;
1246 ATH_MSG_DEBUG("Eff Redefined (to correct for deadfrac): FracDeadStripEta/Phi "
1247 << " Eta/PhiPanelEfficiency " << EtaPanelEfficiency << "/" << PhiPanelEfficiency << " gapEff "
1248 << GapEfficiency);
1249 }
1250 }
1251
1252 // values from COOLDB (eventually overwritten later)
1253 PhiAndEtaEff = float(EtaPanelEfficiency + PhiPanelEfficiency - GapEfficiency);
1254 if (PhiAndEtaEff < 0.) PhiAndEtaEff = 0.;
1255 OnlyEtaEff = float(EtaPanelEfficiency - PhiAndEtaEff);
1256 if (OnlyEtaEff < 0.) OnlyEtaEff = 0.;
1257 OnlyPhiEff = float(PhiPanelEfficiency - PhiAndEtaEff);
1258 if (OnlyPhiEff < 0.) OnlyPhiEff = 0.;
1259
1260 // special patch to be true only when m_Efficiency_fromCOOL=true and /RPC/DQMF/ELEMENT_STATUS tag is
1261 // RPCDQMFElementStatus_2012_Jaunuary_26
1262 bool applySpecialPatch = false;
1263 if (m_EfficiencyPatchForBMShighEta && m_Efficiency_fromCOOL) {
1264 if (m_idHelper->stationName(IdEta) == 3)
1265 {
1266 if (std::abs(m_idHelper->stationEta(IdEta)) == 6 && m_idHelper->doubletR(IdEta) == 1 &&
1267 m_idHelper->doubletZ(IdEta) == 2 && m_idHelper->doubletPhi(IdEta) == 1) {
1268 applySpecialPatch = true;
1269 ATH_MSG_WARNING(
1270 "Applying special patch for BMS at |eta|=6 lowPt plane -dbbZ=2 and dbPhi=1 ... will use default eff. for Id "
1271 << m_idHelper->show_to_string(IdEta));
1272 ATH_MSG_WARNING(
1273 "Applying special patch: THIS HAS TO BE DONE IF /RPC/DQMF/ELEMENT_STATUS tag is "
1274 "RPCDQMFElementStatus_2012_Jaunuary_2");
1275 }
1276 }
1277 }
1278
1279 // if projected tracks number too low or inconsistent values get efficiencies from joboption and overwrite previous values
1280 if (applySpecialPatch || RPC_ProjectedTracksEta < m_CutProjectedTracks || RPC_ProjectedTracksEta > 10000000 ||
1281 EtaPanelEfficiency > 1 || EtaPanelEfficiency < 0 || PhiPanelEfficiency > 1 || PhiPanelEfficiency < 0 || GapEfficiency > 1 ||
1282 GapEfficiency < 0) {
1283 if (index > m_PhiAndEtaEff_A.size() || index > m_OnlyEtaEff_A.size() || index > m_OnlyPhiEff_A.size()) {
1284 THROW_EXCEPTION("Index out of array in Detection Efficiency SideA COOLDB" << index << " stationName = " << stationName);
1285 }
1286 if (RPC_ProjectedTracksEta < m_CutProjectedTracks)
1287 ATH_MSG_DEBUG("# of proj tracks = " << RPC_ProjectedTracksEta << " < cut = " << m_CutProjectedTracks
1288 << " resetting eff. from cool with default(python) values ");
1289
1290 PhiAndEtaEff = m_PhiAndEtaEff_A[index];
1291 OnlyEtaEff = m_OnlyEtaEff_A[index];
1292 OnlyPhiEff = m_OnlyPhiEff_A[index];
1293
1294 if (stationEta < 0) {
1295 if (index > m_PhiAndEtaEff_C.size() || index > m_OnlyEtaEff_C.size() || index > m_OnlyPhiEff_C.size()) {
1296 THROW_EXCEPTION("Index out of array in Detection Efficiency SideC COOLDB" << index << " stationName = " << stationName);
1297 }
1298 PhiAndEtaEff = m_PhiAndEtaEff_C[index];
1299 OnlyEtaEff = m_OnlyEtaEff_C[index];
1300 OnlyPhiEff = m_OnlyPhiEff_C[index];
1301 }
1302
1303 // if (m_applyEffThreshold) {
1304 // gabriele Set efficiency from dead strip fraction instead of nominal value
1305 float effgap = PhiAndEtaEff + OnlyEtaEff + OnlyPhiEff;
1306 float s_EtaPanelEfficiency = 1. - FracDeadStripEta;
1307 float s_PhiPanelEfficiency = 1. - FracDeadStripPhi;
1308 float s_PhiAndEtaEff = s_EtaPanelEfficiency * s_PhiPanelEfficiency / effgap;
1309 if (s_PhiAndEtaEff < PhiAndEtaEff) PhiAndEtaEff = s_PhiAndEtaEff;
1310 float s_OnlyEtaEff = s_EtaPanelEfficiency - PhiAndEtaEff;
1311 float s_OnlyPhiEff = s_PhiPanelEfficiency - PhiAndEtaEff;
1312
1313 if (s_OnlyEtaEff < OnlyEtaEff) OnlyEtaEff = s_OnlyEtaEff;
1314 if (s_OnlyPhiEff < OnlyPhiEff) OnlyPhiEff = s_OnlyPhiEff;
1315 // }
1316 }
1317
1318 float VolEff = PhiAndEtaEff + OnlyEtaEff + OnlyPhiEff;
1319 if (VolEff > maxGeomEff) {
1320 PhiAndEtaEff = (PhiAndEtaEff / VolEff) * maxGeomEff;
1321 OnlyEtaEff = (OnlyEtaEff / VolEff) * maxGeomEff;
1322 OnlyPhiEff = (OnlyPhiEff / VolEff) * maxGeomEff;
1323 }
1324
1325 } // End eff from COOL
1326
1327 // Efficiency correction factor for fractional-charged particles(added by Quanyin Li: quli@cern.ch)
1328 // link to truth particles and calculate the charge and betagamma
1329 HepMC::ConstGenParticlePtr genparticle = trkParticle.cptr();
1330 if (genparticle) {
1331 // only apply efficiency correction to fractional-charged particles based on pdgId betagamma
1332 if (MC::isGenericMultichargedParticle(genparticle)) {
1333 const double eff_sf = FCPEfficiency(genparticle);
1334 // Apply scale factor to the 3 Eff.
1335 PhiAndEtaEff = PhiAndEtaEff * eff_sf;
1336 OnlyEtaEff = OnlyEtaEff * eff_sf;
1337 OnlyPhiEff = OnlyPhiEff * eff_sf;
1338 }
1339 }
1340
1341 float I0 = PhiAndEtaEff;
1342 float I1 = PhiAndEtaEff + OnlyEtaEff;
1343 float ITot = PhiAndEtaEff + OnlyEtaEff + OnlyPhiEff;
1344
1345 float GapEff = ITot ;
1346 float PhiEff = PhiAndEtaEff + OnlyPhiEff;
1347 float EtaEff = PhiAndEtaEff + OnlyEtaEff;
1348
1349 ATH_MSG_DEBUG("DetectionEfficiency: Final Efficiency Values applied for "
1350 << m_idHelper->show_to_string(IdEta) << " are " << PhiAndEtaEff << "=PhiAndEtaEff " << OnlyEtaEff
1351 << "=OnlyEtaEff " << OnlyPhiEff << "=OnlyPhiEff " << GapEff << "=GapEff " << EtaEff << "=EtaEff " << PhiEff
1352 << "=PhiEff ");
1353
1354 float rndmEff = CLHEP::RandFlat::shoot(rndmEngine, 1);
1355
1356 if (rndmEff < I0) {
1357 phiStripOn = true;
1358 etaStripOn = true;
1359 } else if ((I0 <= rndmEff) && (rndmEff < I1)) {
1360 phiStripOn = false;
1361 etaStripOn = true;
1362 } else if ((I1 <= rndmEff) && (rndmEff <= ITot)) {
1363 phiStripOn = true;
1364 etaStripOn = false;
1365 } else {
1366 phiStripOn = false;
1367 etaStripOn = false;
1368 }
1369
1370 return std::make_pair(etaStripOn, phiStripOn);
1371}
1372
1373//--------------------------------------------
1374int RpcDigitizationTool::determineClusterSize(const EventContext& ctx,
1375 const Identifier& idRpcStrip,
1376 double xstripnorm,
1377 CLHEP::HepRandomEngine* rndmEngine) const {
1378 ATH_MSG_DEBUG("RpcDigitizationTool::in determineClusterSize");
1379
1380 ATH_MSG_DEBUG("Digit Id = " << m_idHelper->show_to_string(idRpcStrip));
1381
1382 int ClusterSize = 1;
1383
1384 double FracClusterSize1{1.}, FracClusterSize2{0.}, MeanClusterSize{1.},
1385 FracClusterSizeTail{0.}, MeanClusterSizeTail{1.},
1386 FracClusterSize2norm{0.};
1387
1388 // 2=BML,3=BMS,4=BOL,5=BOS,8=BMF,9=BOF,10=BOG
1389 int stationName = m_idHelper->stationName(idRpcStrip);
1390 int stationEta = m_idHelper->stationEta(idRpcStrip);
1391 int measuresPhi = m_idHelper->measuresPhi(idRpcStrip);
1392
1393 unsigned int index = stationName - 2;
1394 // BML and BMS, BOL and BOS come first (stationName= 2 and 3, 4 and 5 -> index 0-3)
1395 if (stationName > 5 && stationName < 50) index = index - 2;
1396 // BMF, BOF and BOG are 8,9,10 => must be 4,5 and 6
1397 else if (stationName > 50)
1398 index = index - 44;
1399 // BME and BOE 53 and 54 are at indices 7 and 8
1400
1401 if (!m_ClusterSize_fromCOOL && stationName >= 2) {
1402 index += m_FracClusterSize1_A.size() / 2 * measuresPhi;
1403 if (index >= m_FracClusterSize1_A.size() ||
1404 index >= m_FracClusterSize2_A.size() ||
1405 index >= m_FracClusterSizeTail_A.size() ||
1406 index >= m_MeanClusterSizeTail_A.size()) {
1407 ATH_MSG_ERROR("Index out of array in determineClusterSize SideA " << index << " statName " << stationName);
1408 return 1;
1409 }
1410 FracClusterSize1 = m_FracClusterSize1_A[index];
1411 FracClusterSize2 = m_FracClusterSize2_A[index];
1412 FracClusterSizeTail = m_FracClusterSizeTail_A[index];
1413 MeanClusterSizeTail = m_MeanClusterSizeTail_A[index];
1414
1415 if (stationEta < 0) {
1416 index += m_FracClusterSize1_C.size() / 2 * measuresPhi - m_FracClusterSize1_A.size() / 2 * measuresPhi;
1417 if (index >= m_FracClusterSize1_C.size() ||
1418 index >= m_FracClusterSize2_C.size() ||
1419 index >= m_FracClusterSizeTail_C.size() ||
1420 index >= m_MeanClusterSizeTail_C.size()) {
1421 ATH_MSG_ERROR("Index out of array in determineClusterSize SideC " << index << " statName " << stationName);
1422 return 1;
1423 }
1424 FracClusterSize1 = m_FracClusterSize1_C[index];
1425 FracClusterSize2 = m_FracClusterSize2_C[index];
1426 FracClusterSizeTail = m_FracClusterSizeTail_C[index];
1427 MeanClusterSizeTail = m_MeanClusterSizeTail_C[index];
1428 }
1429 } else if (stationName < 2 && (!m_ClusterSize_fromCOOL || !m_ClusterSize_BIS78_fromCOOL)) { // BIS78
1430 FracClusterSize1 = m_FracClusterSize1_BIS78;
1431 FracClusterSize2 = m_FracClusterSize2_BIS78;
1432 FracClusterSizeTail = m_FracClusterSizeTail_BIS78;
1433 MeanClusterSizeTail = m_MeanClusterSizeTail_BIS78;
1434 } else { // Cluster size from COOL
1435 const RpcCondDbData* readCdo{nullptr};
1436 retrieveCondData(ctx, m_readKey, readCdo).ignore();
1437
1438 Identifier Id = m_idHelper->panelID(idRpcStrip);
1439
1440 int RPC_ProjectedTracks = readCdo->getProjectedTrack(Id).value_or(0);
1441 FracClusterSize1 = readCdo->getFracClusterSize1(Id).value_or(1.);
1442 FracClusterSize2 = readCdo->getFracClusterSize2(Id).value_or(0.);
1443 MeanClusterSize = readCdo->getMeanClusterSize(Id).value_or(1.);
1444
1445
1446 ATH_MSG_DEBUG("FracClusterSize1 and 2 " << FracClusterSize1 << " " << FracClusterSize2);
1447
1448 FracClusterSizeTail = 1. - FracClusterSize1 - FracClusterSize2;
1449
1450 MeanClusterSizeTail = MeanClusterSize - FracClusterSize1 - 2 * FracClusterSize2;
1451
1452 ATH_MSG_DEBUG("MeanClusterSizeTail and FracClusterSizeTail " << MeanClusterSizeTail << " " << FracClusterSizeTail);
1453
1454 // if clustersize have anomalous values set to the average cluster size from joboption
1455 if (RPC_ProjectedTracks < m_CutProjectedTracks || RPC_ProjectedTracks > 10000000 || MeanClusterSize > m_CutMaxClusterSize ||
1456 MeanClusterSize <= 1 || FracClusterSizeTail < 0 || FracClusterSize1 < 0 || FracClusterSize2 < 0 || FracClusterSizeTail > 1 ||
1457 FracClusterSize1 > 1 || FracClusterSize2 > 1) {
1458 if (stationName >= 2) {
1459 index += m_FracClusterSize1_A.size() / 2 * measuresPhi;
1460 if (index >= m_FracClusterSize1_A.size() ||
1461 index >= m_FracClusterSize2_A.size() ||
1462 index >= m_FracClusterSizeTail_A.size() ||
1463 index >= m_MeanClusterSizeTail_A.size()) {
1464 ATH_MSG_ERROR("Index out of array in determineClusterSize SideA " << index << " statName " << stationName);
1465 return 1;
1466 }
1467 FracClusterSize1 = m_FracClusterSize1_A[index];
1468 FracClusterSize2 = m_FracClusterSize2_A[index];
1469 FracClusterSizeTail = m_FracClusterSizeTail_A[index];
1470 MeanClusterSizeTail = m_MeanClusterSizeTail_A[index];
1471
1472 if (stationEta < 0) {
1473 index += m_FracClusterSize1_C.size() / 2 * measuresPhi - m_FracClusterSize1_A.size() / 2 * measuresPhi;
1474 if (index > m_FracClusterSize1_C.size() || index > m_FracClusterSize2_C.size() ||
1475 index > m_FracClusterSizeTail_C.size() || index > m_MeanClusterSizeTail_C.size()) {
1476 ATH_MSG_ERROR("Index out of array in determineClusterSize SideC " << index << " statName " << stationName);
1477 return 1;
1478 }
1479
1480 FracClusterSize1 = m_FracClusterSize1_C[index];
1481 FracClusterSize2 = m_FracClusterSize2_C[index];
1482 FracClusterSizeTail = m_FracClusterSizeTail_C[index];
1483 MeanClusterSizeTail = m_MeanClusterSizeTail_C[index];
1484 }
1485 } else {
1486 FracClusterSize1 = m_FracClusterSize1_BIS78;
1487 FracClusterSize2 = m_FracClusterSize2_BIS78;
1488 FracClusterSizeTail = m_FracClusterSizeTail_BIS78;
1489 MeanClusterSizeTail = m_MeanClusterSizeTail_BIS78;
1490 }
1491 }
1492 }
1493 FracClusterSize1 = std::min(FracClusterSize1, 1.);
1494 FracClusterSize2 = std::min(FracClusterSize2, 1.);
1495 FracClusterSizeTail = std::min(FracClusterSizeTail, 1.);
1496 float FracTot = FracClusterSize1 + FracClusterSize2 + FracClusterSizeTail;
1497 if (FracTot != 1. && FracTot > 0) {
1498 FracClusterSize1 = FracClusterSize1 / FracTot;
1499 FracClusterSize2 = FracClusterSize2 / FracTot;
1500 FracClusterSizeTail = FracClusterSizeTail / FracTot;
1501 }
1502 if (MeanClusterSizeTail < 0 || MeanClusterSizeTail > 10) MeanClusterSizeTail = 1;
1503
1504 ATH_MSG_VERBOSE("ClusterSize Final " << FracClusterSize1 << " FracClusterSize1 " << FracClusterSize2 << " FracClusterSize2 "
1505 << FracClusterSizeTail << " " << FracClusterSizeTail << " MeanClusterSizeTail "
1506 << MeanClusterSizeTail);
1507
1508 float FracClusterSize1plus2 = FracClusterSize1 + FracClusterSize2;
1509 float ITot = FracClusterSize1 + FracClusterSize2 + FracClusterSizeTail;
1510
1511 if (FracClusterSize1plus2 != 0) {
1512 // FracClusterSize1norm = FracClusterSize1 / FracClusterSize1plus2 ; // not used
1513 FracClusterSize2norm = FracClusterSize2 / FracClusterSize1plus2;
1514 }
1515
1516 float rndmCS = CLHEP::RandFlat::shoot(rndmEngine, ITot);
1517
1518 if (stationName >= 2) { // Legacy RPCs
1519 // Expanded CS2 of 1.3 to match average CS1 and CS2 (to be investigate)
1520 if (rndmCS < FracClusterSize1plus2) {
1521 // deterministic assignment of CS 1 or 2
1522 if (xstripnorm <= FracClusterSize2norm / 2. * 1.3) {
1523 ClusterSize = -2;
1524 } else if ((1.0 - FracClusterSize2norm / 2. * 1.3) <= xstripnorm) {
1525 ClusterSize = 2;
1526 } else {
1527 ClusterSize = 1;
1528 }
1529 if (m_ClusterSize1_2uncorr) {
1530 float rndmCS1_2 = CLHEP::RandFlat::shoot(rndmEngine, 1);
1531 ClusterSize = 1 + (rndmCS1_2 < FracClusterSize2norm);
1532 }
1533
1534 } else if ((FracClusterSize1plus2 <= rndmCS) && (rndmCS <= ITot)) {
1535 ClusterSize = m_FirstClusterSizeInTail;
1536 ClusterSize += int(CLHEP::RandExponential::shoot(rndmEngine, MeanClusterSizeTail));
1537 float rndmLR = CLHEP::RandFlat::shoot(rndmEngine, 1.0);
1538 if (rndmLR > 0.5) ClusterSize = -ClusterSize;
1539 } else {
1540 ClusterSize = 1;
1541 }
1542
1543 } else { // NRPCs
1544 if (rndmCS < FracClusterSize1) {
1545 ClusterSize = 1;
1546 } else if (rndmCS < FracClusterSize1 + FracClusterSize2) {
1547 ClusterSize = 2;
1548 } else {
1549 ClusterSize = int(CLHEP::RandExponential::shoot(rndmEngine, MeanClusterSizeTail));
1550 }
1551 ClusterSize = std::max(ClusterSize, 1);
1552 if (ClusterSize > 1) {
1553 float rndmLR = CLHEP::RandFlat::shoot(rndmEngine, 1.0);
1554 if (rndmLR > 0.5) ClusterSize = -ClusterSize;
1555 }
1556 }
1557
1558 // negative CS correspond to left asymmetric cluster with respect to nstrip
1559 return ClusterSize;
1560}
1561double RpcDigitizationTool::FCPEfficiency(const HepMC::ConstGenParticlePtr& genParticle) const {
1562 double qcharge = 1.;
1563 const int particlePdgId = genParticle->pdg_id();
1564 // charge calculation
1565 qcharge = (static_cast<double>((std::abs(particlePdgId) / 1000) % 100)) / (static_cast<double>((std::abs(particlePdgId) / 10) % 100));
1566 qcharge = ((static_cast<double>((static_cast<int>(qcharge * 100))))) / 100;
1567 if (particlePdgId < 0.0) qcharge = -qcharge;
1568 // BetaGamma calculation
1569 const double QPx = genParticle->momentum().px();
1570 const double QPy = genParticle->momentum().py();
1571 const double QPz = genParticle->momentum().pz();
1572 const double QE = genParticle->momentum().e();
1573 const double QM2 = std::pow(QE, 2) - std::pow(QPx, 2) - std::pow(QPy, 2) - std::pow(QPz, 2);
1574 const double QP = std::hypot(QPx, QPy, QPz);
1575 const double QM = QM2 >=0 ? std::sqrt(QM2) : -1.;
1576
1577 const double qbetagamma = QM > 0. ? QP / QM : -1.;
1578
1579 // find the i in the array
1580 int i_e = -1;
1581 for (int i = 0; i < 12; i++) {
1582 if (Charge[i] == std::abs(qcharge)) {
1583 i_e = i;
1584 break;
1585 }
1586 }
1587 int i_v = -99, j_v = 99;
1588 if (qbetagamma != -1) {
1589 for (int i = 0; i < 15; i++) {
1590 if (Velocity[i] <= qbetagamma) { i_v = i; }
1591 }
1592 for (int i = 14; i >= 0; i--) {
1593 if (Velocity[i] >= qbetagamma) { j_v = i; }
1594 }
1595 }
1596 // calculate the efficiency according to charge and velocity. Using linear function to calculate efficiency of a specific velocity
1597 // between velocity1 and velocity2
1598 double eff_fcp = 1.0, eff_muon = 1.0;
1599 if (i_e >= 0 && i_e <= 11) {
1600 if (validIndex(j_v, N_Velocity) && validIndex(i_v, N_Velocity) && (j_v - i_v) == 1) {
1601 const double delta_v = Velocity[i_v] - Velocity[j_v];
1602 eff_fcp = (Eff_garfield[i_e][i_v] - Eff_garfield[i_e][j_v]) / delta_v * qbetagamma +
1603 (Eff_garfield[i_e][j_v] * Velocity[i_v] - Eff_garfield[i_e][i_v] * Velocity[j_v]) / delta_v;
1604 eff_muon = (Eff_garfield[11][i_v] - Eff_garfield[11][j_v]) / delta_v * qbetagamma +
1605 (Eff_garfield[11][j_v] * Velocity[i_v] - Eff_garfield[11][i_v] * Velocity[j_v]) / delta_v;
1606 } else if (i_v == 14 && j_v == 99) {
1607 eff_fcp = Eff_garfield[i_e][14];
1608 eff_muon = Eff_garfield[11][14];
1609 } else if (i_v == -99 && j_v == 0) {
1610 eff_fcp = Eff_garfield[i_e][0];
1611 eff_muon = Eff_garfield[11][0];
1612 } else {
1613 ATH_MSG_WARNING("Wrong particle with unknown velocity! Scale factor is set to be 1.");
1614 }
1615 } else {
1616 ATH_MSG_WARNING("Wrong particle with unknown charge! Scale factor is set to be 1.");
1617 }
1618 // A scale factor is calculated by efficiency of fcp / efficiency of muon(charge==1.0
1619 const double eff_SF = eff_fcp / eff_muon;
1620 return eff_SF;
1621}
1622
1623double RpcDigitizationTool::timeOverThreshold(CLHEP::HepRandomEngine* rndmEngine) {
1624 //mn Time-over-threshold modeled as a narrow and a wide gaussian
1625 //mn based on the fit documented in https://its.cern.ch/jira/browse/ATLASRECTS-7820
1626 constexpr double tot_mean_narrow = 16.;
1627 constexpr double tot_sigma_narrow = 2.;
1628 constexpr double tot_mean_wide = 15.;
1629 constexpr double tot_sigma_wide = 4.5;
1630
1631 double thetot = 0.;
1632
1633 if (CLHEP::RandFlat::shoot(rndmEngine)<0.75) {
1634 thetot = CLHEP::RandGaussZiggurat::shoot(rndmEngine, tot_mean_narrow, tot_sigma_narrow);
1635 } else {
1636 thetot = CLHEP::RandGaussZiggurat::shoot(rndmEngine, tot_mean_wide, tot_sigma_wide);
1637 }
1638
1639 return (thetot > 0.) ? thetot : 0.;
1640}
hitTime
float hitTime(const AFP_SIDSimHit &hit)
Definition AFP_SIDSimHit.h:39
make_array
constexpr std::array< T, N > make_array(const T &def_val)
Helper function to initialize in-place arrays with non-zero values.
Definition ArrayHelper.h:10
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition AthMsgStreamMacros.h:33
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition AthMsgStreamMacros.h:34
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition AthMsgStreamMacros.h:31
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition AthMsgStreamMacros.h:28
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
isValid
bool isValid(const T &p)
Av: we implement here an ATLAS-sepcific convention: all particles which are 99xxxxx are fine.
Definition AtlasPID.h:878
HepMCHelpers.h
ATLAS-specific HepMC functions.
SubEventIterator
std::vector< xAOD::EventInfo::SubEvent >::const_iterator SubEventIterator
Definition IPileUpTool.h:22
IRDBAccessSvc.h
Definition of the abstract IRDBAccessSvc interface.
IRDBRecordset_ptr
std::shared_ptr< IRDBRecordset > IRDBRecordset_ptr
Definition IRDBAccessSvc.h:27
IRDBRecord.h
Definition of the abstract IRDBRecord interface.
IRDBRecordset.h
Definition of the abstract IRDBRecordset interface.
sc
static Double_t sc
Definition LArPhysWaveHECTool.cxx:37
RPCSimHitCollection
AtlasHitsVector< RPCSimHit > RPCSimHitCollection
Definition RPCSimHitCollection.h:15
Unknown
@ Unknown
Definition TruthClasses.h:9
x
#define x
ATHRNG::RNGWrapper
A wrapper class for event-slot-local random engines.
Definition RNGWrapper.h:56
ATHRNG::RNGWrapper::setSeed
void setSeed(const std::string &algName, const EventContext &ctx)
Set the random seed using a string (e.g.
Definition RNGWrapper.h:169
ATHRNG::RNGWrapper::getEngine
CLHEP::HepRandomEngine * getEngine(const EventContext &ctx) const
Retrieve the random engine corresponding to the provided EventContext.
Definition RNGWrapper.h:134
AtlasHitsVector::size
size_type size() const
Definition AtlasHitsVector.h:142
DataVector::get
const T * get(size_type n) const
Access an element, as an rvalue.
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
IRDBRecordset::size
virtual unsigned int size() const =0
IdContext
This class saves the "context" of an expanded identifier (ExpandedIdentifier) for compact or hash ver...
Definition IdContext.h:26
IdContext::begin_index
size_type begin_index() const
Definition IdContext.h:45
IdContext::end_index
size_type end_index() const
Definition IdContext.h:46
Identifier32::get_compact
value_type get_compact() const
Get the compact id.
Definition Identifier32.h:47
IdentifierHash
This is a "hash" representation of an Identifier.
Definition IdentifierHash.h:25
Identifier::show
void show(std::ostream &out=std::cout) const
Print out in hex form.
Definition Identifier.cxx:30
Identifier::get_identifier32
Identifier32 get_identifier32() const
Get the 32-bit version Identifier, will be invalid if >32 bits needed.
MuonGM::MuonClusterReadoutElement::transform
virtual const Amg::Transform3D & transform() const override
Return local to global transform.
Definition MuonClusterReadoutElement.h:124
MuonGM::MuonDetectorManager
The MuonDetectorManager stores the transient representation of the Muon Spectrometer geometry and pro...
Definition MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonDetectorManager.h:51
MuonGM::MuonDetectorManager::getRpcReadoutElement
const RpcReadoutElement * getRpcReadoutElement(const Identifier &id) const
access via extended identifier (requires unpacking)
Definition MuonDetDescr/MuonReadoutGeometry/src/MuonDetectorManager.cxx:180
MuonGM::RpcReadoutElement
An RpcReadoutElement corresponds to a single RPC module; therefore typicaly a barrel muon station con...
Definition MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/RpcReadoutElement.h:55
MuonGM::RpcReadoutElement::nGasGapPerLay
int nGasGapPerLay() const
returns the number of gasgaps
MuonGM::RpcReadoutElement::NphiStripPanels
int NphiStripPanels() const
returns the number of phi strip panels (1 or 2)
MuonGM::RpcReadoutElement::stripNumber
virtual int stripNumber(const Amg::Vector2D &pos, const Identifier &id) const override final
strip number corresponding to local position.
MuonGM::RpcReadoutElement::localGasGapPos
Amg::Vector3D localGasGapPos(const Identifier &id) const
Returns the position of the gasGap w.r.t. rest frame of the chamber.
Definition MuonDetDescr/MuonReadoutGeometry/src/RpcReadoutElement.cxx:193
MuonGM::RpcReadoutElement::localToGlobalTransf
Amg::Transform3D localToGlobalTransf(const Identifier &id) const
Definition MuonDetDescr/MuonReadoutGeometry/src/RpcReadoutElement.cxx:211
MuonGM::RpcReadoutElement::Nstrips
int Nstrips(bool measphi) const
returns the number of strips for the phi or eta plane
MuonGM::RpcReadoutElement::numberOfLayers
virtual int numberOfLayers(bool measphi=true) const override final
number of layers in phi/eta projection, same for eta/phi planes
MuonGM::RpcReadoutElement::StripPitch
double StripPitch(bool measphi) const
returns the strip pitch for the phi or eta plane
MuonSimData::Deposit
std::pair< HepMcParticleLink, MuonMCData > Deposit
Definition MuonSimData.h:66
PileUpToolBase::m_vetoPileUpTruthLinks
Gaudi::Property< int > m_vetoPileUpTruthLinks
Definition PileUpToolBase.h:58
PileUpToolBase::PileUpToolBase
PileUpToolBase(const std::string &type, const std::string &name, const IInterface *parent)
Definition PileUpToolBase.cxx:7
RPCSimHit::particleEncoding
int particleEncoding() const
Definition RPCSimHit.h:61
RPCSimHit::stepLength
double stepLength() const
Definition RPCSimHit.h:60
RPCSimHit::postLocalPosition
const Amg::Vector3D & postLocalPosition() const
Definition RPCSimHit.h:57
RPCSimHit::RPCid
int RPCid() const
Definition RPCSimHit.h:63
RPCSimHit::localPosition
const Amg::Vector3D & localPosition() const
Definition RPCSimHit.h:55
RPCSimHit::kineticEnergy
double kineticEnergy() const
Definition RPCSimHit.h:62
RPCSimHit::globalTime
double globalTime() const
Definition RPCSimHit.h:54
RPCSimHit::energyDeposit
double energyDeposit() const
Definition RPCSimHit.h:58
RpcCondDbData::getFracClusterSize2
std::optional< double > getFracClusterSize2(const Identifier &) const
Definition RpcCondDbData.cxx:44
RpcCondDbData::getProjectedTrack
std::optional< int > getProjectedTrack(const Identifier &) const
Definition RpcCondDbData.cxx:60
RpcCondDbData::getFracDeadStrip
std::optional< double > getFracDeadStrip(const Identifier &) const
Definition RpcCondDbData.cxx:48
RpcCondDbData::getFracClusterSize1
std::optional< double > getFracClusterSize1(const Identifier &) const
Definition RpcCondDbData.cxx:40
RpcCondDbData::getGapEfficiency
std::optional< double > getGapEfficiency(const Identifier &) const
Definition RpcCondDbData.cxx:52
RpcCondDbData::getEfficiency
std::optional< double > getEfficiency(const Identifier &) const
Definition RpcCondDbData.cxx:36
RpcCondDbData::getMeanClusterSize
std::optional< double > getMeanClusterSize(const Identifier &) const
Definition RpcCondDbData.cxx:56
RpcDigitizationTool::detectionEfficiency
std::pair< bool, bool > detectionEfficiency(const EventContext &ctx, const Identifier &ideta, const Identifier &idphi, CLHEP::HepRandomEngine *rndmEngine, const HepMcParticleLink &trkParticle) const
Evaluate detection efficiency.
Definition RpcDigitizationTool.cxx:1088
RpcDigitizationTool::m_MeanClusterSizeTail_BIS78
Gaudi::Property< float > m_MeanClusterSizeTail_BIS78
Definition RpcDigitizationTool.h:204
RpcDigitizationTool::m_CorrJitter_BIS78
Gaudi::Property< double > m_CorrJitter_BIS78
Definition RpcDigitizationTool.h:130
RpcDigitizationTool::Collections_t
std::vector< std::unique_ptr< RpcDigitCollection > > Collections_t
Definition RpcDigitizationTool.h:94
RpcDigitizationTool::outsideWindow
bool outsideWindow(double time) const
Definition RpcDigitizationTool.h:257
RpcDigitizationTool::m_muonHelper
const RpcHitIdHelper * m_muonHelper
Definition RpcDigitizationTool.h:161
RpcDigitizationTool::determineClusterSize
int determineClusterSize(const EventContext &ctx, const Identifier &id, double xstripnorm, CLHEP::HepRandomEngine *rndmEngine) const
Definition RpcDigitizationTool.cxx:1374
RpcDigitizationTool::mergeEvent
virtual StatusCode mergeEvent(const EventContext &ctx) override final
When being run from PileUpToolsAlgs, this method is called at the end of the subevts loop.
Definition RpcDigitizationTool.cxx:405
RpcDigitizationTool::initialize
virtual StatusCode initialize() override final
Initialize.
Definition RpcDigitizationTool.cxx:96
RpcDigitizationTool::PackMCTruth
long long int PackMCTruth(float proptime, float tof, float posx, float posz) const
Definition RpcDigitizationTool.cxx:1008
RpcDigitizationTool::m_turnON_efficiency
Gaudi::Property< bool > m_turnON_efficiency
Definition RpcDigitizationTool.h:175
RpcDigitizationTool::FCPEfficiency
double FCPEfficiency(const HepMC::ConstGenParticlePtr &genParticle) const
Definition RpcDigitizationTool.cxx:1561
RpcDigitizationTool::m_EfficiencyPatchForBMShighEta
Gaudi::Property< bool > m_EfficiencyPatchForBMShighEta
Definition RpcDigitizationTool.h:231
RpcDigitizationTool::m_patch_for_rpc_time
Gaudi::Property< bool > m_patch_for_rpc_time
Definition RpcDigitizationTool.h:167
RpcDigitizationTool::physicalClusterSize
std::array< int, 3 > physicalClusterSize(const EventContext &ctx, const MuonGM::RpcReadoutElement *reEle, const Identifier &id, const Amg::Vector3D &posAtCentre, CLHEP::HepRandomEngine *rndmEngine) const
Cluster simulation: first step.
Definition RpcDigitizationTool.cxx:922
RpcDigitizationTool::RpcDigitizationTool
RpcDigitizationTool(const std::string &type, const std::string &name, const IInterface *pIID)
Definition RpcDigitizationTool.cxx:91
RpcDigitizationTool::m_FirstClusterSizeInTail
Gaudi::Property< int > m_FirstClusterSizeInTail
Definition RpcDigitizationTool.h:178
RpcDigitizationTool::m_FracClusterSize2_A
Gaudi::Property< std::vector< double > > m_FracClusterSize2_A
Definition RpcDigitizationTool.h:192
RpcDigitizationTool::m_rndmSvc
ServiceHandle< IAthRNGSvc > m_rndmSvc
Definition RpcDigitizationTool.h:224
RpcDigitizationTool::m_timeWindowLowerOffset
Gaudi::Property< double > m_timeWindowLowerOffset
Definition RpcDigitizationTool.h:140
RpcDigitizationTool::m_FracClusterSize2_BIS78
Gaudi::Property< float > m_FracClusterSize2_BIS78
Definition RpcDigitizationTool.h:202
RpcDigitizationTool::fromSimHitToLayer
Amg::Transform3D fromSimHitToLayer(const MuonGM::RpcReadoutElement *readOutEle, const Identifier &layerId) const
Returns the position of the hit expressed in the gasGap coordinate system.
Definition RpcDigitizationTool.cxx:900
RpcDigitizationTool::m_CutProjectedTracks
Gaudi::Property< int > m_CutProjectedTracks
Definition RpcDigitizationTool.h:246
RpcDigitizationTool::m_timeWindowUpperOffset
Gaudi::Property< double > m_timeWindowUpperOffset
Definition RpcDigitizationTool.h:141
RpcDigitizationTool::m_FracClusterSize1_A
Gaudi::Property< std::vector< double > > m_FracClusterSize1_A
Definition RpcDigitizationTool.h:191
RpcDigitizationTool::m_deadTime
Gaudi::Property< int > m_deadTime
Definition RpcDigitizationTool.h:166
RpcDigitizationTool::m_kill_deadstrips
Gaudi::Property< bool > m_kill_deadstrips
Definition RpcDigitizationTool.h:176
RpcDigitizationTool::m_ClusterSize_fromCOOL
Gaudi::Property< bool > m_ClusterSize_fromCOOL
Definition RpcDigitizationTool.h:234
RpcDigitizationTool::TurnOnStrips
std::array< int, 3 > TurnOnStrips(const MuonGM::RpcReadoutElement *reEle, std::array< int, 3 > &&pcs, const Identifier &id) const
Cluster simulation: second step.
Definition RpcDigitizationTool.cxx:961
RpcDigitizationTool::m_includePileUpTruth
Gaudi::Property< bool > m_includePileUpTruth
Definition RpcDigitizationTool.h:173
RpcDigitizationTool::m_simHitValidKey
SG::WriteHandleKey< RPCSimHitCollection > m_simHitValidKey
Definition RpcDigitizationTool.h:223
RpcDigitizationTool::m_outputDigitCollectionKey
SG::WriteHandleKey< RpcDigitContainer > m_outputDigitCollectionKey
Definition RpcDigitizationTool.h:218
RpcDigitizationTool::initializeRunDependentParameters
StatusCode initializeRunDependentParameters()
Definition RpcDigitizationTool.cxx:176
RpcDigitizationTool::m_ClusterSize_BIS78_fromCOOL
Gaudi::Property< bool > m_ClusterSize_BIS78_fromCOOL
Definition RpcDigitizationTool.h:241
RpcDigitizationTool::m_ClusterSize1_2uncorr
Gaudi::Property< bool > m_ClusterSize1_2uncorr
Definition RpcDigitizationTool.h:235
RpcDigitizationTool::m_ignoreRunDepConfig
Gaudi::Property< bool > m_ignoreRunDepConfig
Definition RpcDigitizationTool.h:238
RpcDigitizationTool::m_mergeSvc
ServiceHandle< PileUpMergeSvc > m_mergeSvc
Definition RpcDigitizationTool.h:213
RpcDigitizationTool::m_MeanClusterSizeTail_A
Gaudi::Property< std::vector< double > > m_MeanClusterSizeTail_A
Definition RpcDigitizationTool.h:194
RpcDigitizationTool::m_RPCHitCollList
std::vector< std::unique_ptr< RPCSimHitCollection > > m_RPCHitCollList
Definition RpcDigitizationTool.h:162
RpcDigitizationTool::m_outputSDO_CollectionKey
SG::WriteHandleKey< MuonSimDataCollection > m_outputSDO_CollectionKey
Definition RpcDigitizationTool.h:220
RpcDigitizationTool::m_readKey
SG::ReadCondHandleKey< RpcCondDbData > m_readKey
Definition RpcDigitizationTool.h:164
RpcDigitizationTool::m_MeanClusterSizeTail_C
Gaudi::Property< std::vector< double > > m_MeanClusterSizeTail_C
Definition RpcDigitizationTool.h:199
RpcDigitizationTool::m_sdo_tmp_map
std::map< Identifier, std::vector< MuonSimData::Deposit > > m_sdo_tmp_map
Definition RpcDigitizationTool.h:165
RpcDigitizationTool::m_OnlyPhiEff_C
Gaudi::Property< std::vector< float > > m_OnlyPhiEff_C
Definition RpcDigitizationTool.h:184
RpcDigitizationTool::m_validationSetup
Gaudi::Property< bool > m_validationSetup
Definition RpcDigitizationTool.h:172
RpcDigitizationTool::m_FracClusterSize1_C
Gaudi::Property< std::vector< double > > m_FracClusterSize1_C
Definition RpcDigitizationTool.h:196
RpcDigitizationTool::m_detMgrKey
SG::ReadCondHandleKey< MuonGM::MuonDetectorManager > m_detMgrKey
Definition RpcDigitizationTool.h:158
RpcDigitizationTool::prepareEvent
virtual StatusCode prepareEvent(const EventContext &ctx, const unsigned int) override final
When being run from PileUpToolsAlgs, this method is called at the start of the subevts loop.
Definition RpcDigitizationTool.cxx:304
RpcDigitizationTool::m_UncorrJitter_BIS78
Gaudi::Property< double > m_UncorrJitter_BIS78
Definition RpcDigitizationTool.h:129
RpcDigitizationTool::processBunchXing
virtual StatusCode processBunchXing(int bunchXing, SubEventIterator bSubEvents, SubEventIterator eSubEvents) override final
When being run from PileUpToolsAlgs, this method is called for each active bunch-crossing to process ...
Definition RpcDigitizationTool.cxx:316
RpcDigitizationTool::PropagationTime
double PropagationTime(const MuonGM::RpcReadoutElement *reEle, const Identifier &id, const Amg::Vector3D &globPos) const
Calculates the propagation time along the strip.
Definition RpcDigitizationTool.cxx:992
RpcDigitizationTool::m_OnlyPhiEff_A
Gaudi::Property< std::vector< float > > m_OnlyPhiEff_A
Definition RpcDigitizationTool.h:181
RpcDigitizationTool::m_FracClusterSizeTail_BIS78
Gaudi::Property< float > m_FracClusterSizeTail_BIS78
Definition RpcDigitizationTool.h:203
RpcDigitizationTool::m_FracClusterSizeTail_A
Gaudi::Property< std::vector< double > > m_FracClusterSizeTail_A
Definition RpcDigitizationTool.h:193
RpcDigitizationTool::m_OnlyEtaEff_BIS78
Gaudi::Property< float > m_OnlyEtaEff_BIS78
Definition RpcDigitizationTool.h:188
RpcDigitizationTool::m_FracClusterSize2_C
Gaudi::Property< std::vector< double > > m_FracClusterSize2_C
Definition RpcDigitizationTool.h:197
RpcDigitizationTool::m_thpcRPC
std::unique_ptr< TimedHitCollection< RPCSimHit > > m_thpcRPC
Definition RpcDigitizationTool.h:163
RpcDigitizationTool::getNextEvent
StatusCode getNextEvent(const EventContext &ctx)
Get next event and extract collection of hit collections:
Definition RpcDigitizationTool.cxx:350
RpcDigitizationTool::m_OnlyEtaEff_C
Gaudi::Property< std::vector< float > > m_OnlyEtaEff_C
Definition RpcDigitizationTool.h:185
RpcDigitizationTool::m_RPC_TimeSchema
Gaudi::Property< std::string > m_RPC_TimeSchema
Definition RpcDigitizationTool.h:226
RpcDigitizationTool::m_idHelper
const RpcIdHelper * m_idHelper
Definition RpcDigitizationTool.h:160
RpcDigitizationTool::UnPackMCTruth
static void UnPackMCTruth(double theWord, float &proptime, float &tof, float &posy, float &posz)
Definition RpcDigitizationTool.cxx:1036
RpcDigitizationTool::timeOverThreshold
static double timeOverThreshold(CLHEP::HepRandomEngine *rndmEngine)
Definition RpcDigitizationTool.cxx:1623
RpcDigitizationTool::m_CutMaxClusterSize
Gaudi::Property< float > m_CutMaxClusterSize
Definition RpcDigitizationTool.h:245
RpcDigitizationTool::m_hitsContainerKey
SG::ReadHandleKey< RPCSimHitCollection > m_hitsContainerKey
Definition RpcDigitizationTool.h:216
RpcDigitizationTool::m_OnlyPhiEff_BIS78
Gaudi::Property< float > m_OnlyPhiEff_BIS78
Definition RpcDigitizationTool.h:189
RpcDigitizationTool::m_turnON_clustersize
Gaudi::Property< bool > m_turnON_clustersize
Definition RpcDigitizationTool.h:177
RpcDigitizationTool::m_PhiAndEtaEff_A
Gaudi::Property< std::vector< float > > m_PhiAndEtaEff_A
Definition RpcDigitizationTool.h:180
RpcDigitizationTool::m_RPCInfoFromDb
Gaudi::Property< bool > m_RPCInfoFromDb
Definition RpcDigitizationTool.h:244
RpcDigitizationTool::m_rpc_time_shift
Gaudi::Property< double > m_rpc_time_shift
Definition RpcDigitizationTool.h:168
RpcDigitizationTool::m_FracClusterSize1_BIS78
Gaudi::Property< float > m_FracClusterSize1_BIS78
Definition RpcDigitizationTool.h:201
RpcDigitizationTool::m_CorrJitter
Gaudi::Property< double > m_CorrJitter
Definition RpcDigitizationTool.h:127
RpcDigitizationTool::m_OnlyEtaEff_A
Gaudi::Property< std::vector< float > > m_OnlyEtaEff_A
Definition RpcDigitizationTool.h:182
RpcDigitizationTool::m_Efficiency_fromCOOL
Gaudi::Property< bool > m_Efficiency_fromCOOL
Definition RpcDigitizationTool.h:230
RpcDigitizationTool::m_UncorrJitter
Gaudi::Property< double > m_UncorrJitter
Calculates the position of the hit wrt to the strip panel this transformation is needed since the imp...
Definition RpcDigitizationTool.h:126
RpcDigitizationTool::m_inputHitCollectionName
std::string m_inputHitCollectionName
Definition RpcDigitizationTool.h:217
RpcDigitizationTool::doDigitization
StatusCode doDigitization(const EventContext &ctx, Collections_t &collections, MuonSimDataCollection *sdoContainer)
Digitization functionality shared with RPC_PileUpTool.
Definition RpcDigitizationTool.cxx:480
RpcDigitizationTool::processAllSubEvents
virtual StatusCode processAllSubEvents(const EventContext &ctx) override final
alternative interface which uses the PileUpMergeSvc to obtain all the required SubEvents.
Definition RpcDigitizationTool.cxx:439
RpcDigitizationTool::retrieveCondData
StatusCode retrieveCondData(const EventContext &ctx, const SG::ReadCondHandleKey< CondType > &key, const CondType *&condPtr) const
Definition RpcDigitizationTool.cxx:285
RpcDigitizationTool::m_onlyUseContainerName
Gaudi::Property< bool > m_onlyUseContainerName
Definition RpcDigitizationTool.h:214
RpcDigitizationTool::m_PhiAndEtaEff_C
Gaudi::Property< std::vector< float > > m_PhiAndEtaEff_C
Definition RpcDigitizationTool.h:183
RpcDigitizationTool::m_Efficiency_BIS78_fromCOOL
Gaudi::Property< bool > m_Efficiency_BIS78_fromCOOL
Definition RpcDigitizationTool.h:240
RpcDigitizationTool::m_FracClusterSizeTail_C
Gaudi::Property< std::vector< double > > m_FracClusterSizeTail_C
Definition RpcDigitizationTool.h:198
RpcDigitizationTool::m_BOG_BOF_DoubletR2_OFF
Gaudi::Property< bool > m_BOG_BOF_DoubletR2_OFF
Definition RpcDigitizationTool.h:237
RpcDigitizationTool::m_muonOnlySDOs
Gaudi::Property< bool > m_muonOnlySDOs
Definition RpcDigitizationTool.h:208
RpcDigitizationTool::m_sdoAreOnlyDigits
Gaudi::Property< bool > m_sdoAreOnlyDigits
Definition RpcDigitizationTool.h:227
RpcDigitizationTool::m_PhiAndEtaEff_BIS78
Gaudi::Property< float > m_PhiAndEtaEff_BIS78
Definition RpcDigitizationTool.h:187
RpcHitIdHelper::GetHelper
static const RpcHitIdHelper * GetHelper(unsigned int nGasGaps=2)
Definition RpcHitIdHelper.cxx:23
SG::ReadCondHandle::cptr
const_pointer_type cptr()
Definition ReadCondHandle.h:67
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
SG::ReadHandle::cptr
const_pointer_type cptr()
Dereference the pointer.
SG::VarHandleBase::store
std::string store() const
Return the name of the store holding the object we are proxying.
Definition StoreGate/src/VarHandleBase.cxx:382
SG::VarHandleBase::name
const std::string & name() const
Return the StoreGate ID for the referenced object.
Definition AthToolSupport/AsgDataHandles/Root/VarHandleBase.cxx:75
SG::WriteHandle::record
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
SG::WriteHandle::ptr
pointer_type ptr()
Dereference the pointer.
TimedHitCollection::const_iterator
TimedVector::const_iterator const_iterator
Definition TimedHitCollection.h:20
TimedHitPtr
a smart pointer to a hit that also provides access to the extended timing info of the host event.
Definition TimedHitPtr.h:18
TimedHitPtr::eventId
unsigned short eventId() const
the index of the component event in PileUpEventInfo.
Definition TimedHitPtr.h:47
simData
constexpr bool simData
Definition constants.h:36
Amg::getTranslate3D
Amg::Transform3D getTranslate3D(const double X, const double Y, const double Z)
: Returns a shift transformation along an arbitrary axis
Definition GeoPrimitivesHelpers.h:289
Amg::intersect
std::optional< double > intersect(const AmgVector(N)&posA, const AmgVector(N)&dirA, const AmgVector(N)&posB, const AmgVector(N)&dirB)
Calculates the point B' along the line B that's closest to a second line A.
Definition GeoPrimitivesHelpers.h:347
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition GeoPrimitivesToStringConverter.h:40
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition GeoPrimitives.h:46
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
Amg::getRotateY3D
Amg::Transform3D getRotateY3D(double angle)
get a rotation transformation around Y-axis
Definition GeoPrimitivesHelpers.h:261
HepMC::ConstGenParticlePtr
const GenParticle * ConstGenParticlePtr
Definition GenParticle.h:38
HepMC::ignoreTruthLink
bool ignoreTruthLink(const T &p, bool vetoPileUp)
Helper function for SDO creation in PileUpTools.
Definition MagicNumbers.h:345
MC::isGenericMultichargedParticle
bool isGenericMultichargedParticle(const T &p)
In addition, there is a need to identify ”Q-ball” and similar very exotic (multi-charged) particles w...
Definition HepMCHelpers.h:685
MuonGM
Ensure that the Athena extensions are properly loaded.
Definition GeoMuonHits.h:27
Run1
USAGE: openCoraCool.exe "COOLONL_SCT/COMP200".
Definition openCoraCool.cxx:57
index
Definition index.py:1
PileUpMergeSvc::TimedList::type
std::list< value_t > type
type of the collection of timed data object
Definition PileUpMergeSvc.h:73
PileUpTimeEventIndex
a struct encapsulating the identifier of a pile-up event
Definition PileUpTimeEventIndex.h:12
PileUpTimeEventIndex::index
index_type index() const
the index of the component event in PileUpEventInfo
Definition PileUpTimeEventIndex.cxx:76
PileUpTimeEventIndex::type
PileUpType type() const
the pileup type - minbias, cavern, beam halo, signal?
Definition PileUpTimeEventIndex.cxx:81
PileUpTimeEventIndex::time
time_type time() const
bunch xing time in ns
Definition PileUpTimeEventIndex.cxx:71
result
type
value_type
Definition EDM_MasterSearch.h:11
THROW_EXCEPTION
#define THROW_EXCEPTION(MESSAGE)
Definition throwExcept.h:10
run
int run(int argc, char *argv[])
Definition ttree2hdf5.cxx:28
Generated on for ATLAS Offline Software by doxygen 1.14.0