10 #include "GaudiKernel/SystemOfUnits.h"
33 #ifdef TRKDETDESCR_MEMUSAGE
39 m_checkForEmptyHits(true),
40 m_mappingVolumeName(
"Atlas"),
41 m_mappingVolume(nullptr),
42 m_inputMaterialStepCollection(
"MaterialStepRecords"),
45 m_useLayerThickness(false),
48 m_mapComposition(false),
49 m_minCompositionFraction(0.005),
50 m_elementTable(nullptr),
51 m_inputEventElementTable(
"ElementTable"),
52 m_accumulatedMaterialXX0(0.),
53 m_accumulatedRhoS(0.),
57 m_layerMaterialScreenOutput(0)
58 #ifdef TRKDETDESCR_MEMUSAGE
94 ATH_CHECK( m_trackingGeometryReadKey.initialize(!m_trackingGeometryReadKey.key().empty()) );
96 ATH_CHECK(m_extrapolationEngine.retrieve());
98 if ( !m_materialMapper.empty() )
101 if ( !m_layerMaterialRecordAnalyser.empty() )
102 ATH_CHECK( m_layerMaterialRecordAnalyser.retrieve() );
104 if ( !m_layerMaterialCreators.empty() )
105 ATH_CHECK( m_layerMaterialCreators.retrieve() );
107 if ( !m_layerMaterialAnalysers.empty() )
108 ATH_CHECK( m_layerMaterialAnalysers.retrieve() );
110 ATH_CHECK( m_inputMaterialStepCollection.initialize() );
111 ATH_CHECK( m_inputEventElementTable.initialize() );
113 return StatusCode::SUCCESS;
122 if (!m_mappingVolume) {
123 StatusCode retrieveCode = handleTrackingGeometry();
124 if (retrieveCode.isFailure()){
125 ATH_MSG_INFO(
"Could not retrieve mapping volume from tracking geometry. Exiting.");
130 ATH_MSG_VERBOSE(
"Mapping volume correctly retrieved from tracking geometry");
137 if (m_mapComposition) {
139 (*m_elementTable) += (*eTableEvent);
144 int associatedSteps = 0;
145 m_accumulatedMaterialXX0 = 0.;
146 m_accumulatedRhoS = 0.;
147 m_layersRecordedPerEvent.clear();
149 if (materialStepCollection.
isValid() && !materialStepCollection->empty()){
152 size_t materialSteps = materialStepCollection->size();
153 ATH_MSG_DEBUG(
"[+] Successfully read "<< materialSteps <<
" geantino steps");
156 double dirx = (*materialStepCollection)[materialSteps-1]->hitX();
157 double diry = (*materialStepCollection)[materialSteps-1]->hitY();
158 double dirz = (*materialStepCollection)[materialSteps-1]->hitZ();
161 double eta = direction.eta();
163 if ( fabs(eta) > m_etaCutOff || (m_etaSide && m_etaSide*eta < 0.) ) {
164 ATH_MSG_VERBOSE(
"[-] Event is outside eta acceptance of " << m_etaCutOff <<
". Skipping it.");
165 return StatusCode::SUCCESS;
178 std::vector< std::pair<const Trk::Layer*, Amg::Vector3D> > layersAndHits;
185 ATH_MSG_VERBOSE(
"[+] Extrapolation to layers did succeed and found " << nLayersHit <<
" layers.");
187 layersAndHits.reserve(nLayersHit);
210 std::pair<const Trk::Layer*, Amg::Vector3D> layerHitPair(
213 << ++ilayer <<
" with index "
217 layersAndHits.push_back(layerHitPair);
226 if (layersAndHits.empty()){
228 return StatusCode::SUCCESS;
232 size_t currentLayer = 0;
236 ATH_MSG_VERBOSE(
"[L] starting from layer " << currentLayer <<
" from layer collection for this step.");
238 double t =
step->steplength();
244 if (!m_mappingVolume || !(m_mappingVolume->inside(
pos)) || fabs(
pos.eta()) > m_etaCutOff ){
250 if (currentLayer < nLayersHit-1) {
252 double currentDistance = (
pos-layersAndHits[currentLayer].second).
mag();
254 for (
size_t testLayer = (currentLayer+1); testLayer < nLayersHit; ++testLayer){
256 double testDistance = (
pos-layersAndHits[testLayer].second).
mag();
257 ATH_MSG_VERBOSE(
"[L] Testing layer " << testLayer <<
" from layer collection for this step.");
259 if ( testDistance < currentDistance ){
261 ATH_MSG_VERBOSE(
"[L] Skipping over to current layer " << testLayer <<
" because " << testDistance <<
" < " << currentDistance);
263 currentLayer = testLayer;
264 currentDistance = testDistance;
272 const Trk::Layer* assignedLayer = layersAndHits[currentLayer].first;
273 Amg::Vector3D assignedPosition = layersAndHits[currentLayer].second;
278 associateHit(*assignedLayer,
pos, assignedPosition,
t,
step->fullMaterial());
282 ATH_MSG_VERBOSE(
"Found " << layersAndHits.size() <<
" intersected layers - while having " << m_layersRecordedPerEvent.size() <<
" recorded ones.");
285 for (
auto& lhp : layersAndHits){
287 if (m_layersRecordedPerEvent.find(lhp.first) == m_layersRecordedPerEvent.end()){
289 auto clIter = m_layerRecords.find(lhp.first);
290 if (clIter != m_layerRecords.end() ){
291 (*clIter).second.associateEmptyHit(lhp.second);
292 ATH_MSG_VERBOSE(
"- to layer with index "<< lhp.first->layerIndex().value() <<
" with empty hit detected.");
294 ATH_MSG_WARNING(
"- no Layer found in the associated map! Should not happen.");
299 if (associatedSteps) {
300 ATH_MSG_VERBOSE(
"There are associated steps, need to call finalizeEvent() & record to the MaterialMapper.");
302 for (
auto& lRecord : m_layerRecords ) {
306 if (assMatHit && !m_materialMapper.empty()) m_materialMapper->recordLayerHit(*assMatHit,
true);
318 return StatusCode::SUCCESS;
337 auto clIter = m_layerRecords.find(
layer);
338 if (clIter != m_layerRecords.end() ){
340 m_layersRecordedPerEvent[
layer] =
true;
342 (*clIter).second.associateGeantinoHit(positionOnLayer, stepl,
mat);
343 ATH_MSG_VERBOSE(
"- associate Geantino Information at intersection [r/z] = " << positionOnLayer.perp() <<
"/"<< positionOnLayer.z() );
345 ATH_MSG_VERBOSE(
"- associate Geantino Information ( s, s/x0 , x0 , l0, a, z, rho ) = "
346 << stepl <<
", "<< stepl/
mat.X0 <<
", "<<
mat.X0 <<
", "<<
mat.L0 <<
", "<<
mat.A <<
", "<<
mat.Z <<
", "<<
mat.rho );
347 m_accumulatedMaterialXX0 += stepl/
mat.X0;
348 m_accumulatedRhoS += stepl*
mat.rho;
349 ATH_MSG_VERBOSE(
"- accumulated material information X/x0 = " << m_accumulatedMaterialXX0);
350 ATH_MSG_VERBOSE(
"- accumulated effective densitity rho*s = " << m_accumulatedRhoS);
353 if (m_mapMaterial && !m_materialMapper.empty()){
355 m_materialMapper->recordMaterialHit(am, positionOnLayer);
360 ATH_MSG_WARNING(
"- associate hit - the layer with index " <<
layer->layerIndex().value() <<
" was not found - should not happen!");
371 if (!propSet)
return;
377 if (confinedLayers) {
385 if (
layer && (*layer).layerIndex().value() ) {
386 ATH_MSG_INFO(
" > LayerIndex: "<< (*layer).layerIndex() );
390 auto curIt = propSet->find((*layer).layerIndex());
391 if (curIt != propSet->end()) {
392 ATH_MSG_INFO(
"LayerMaterial assigned for Layer with index: "<< (*layer).layerIndex() );
394 layer->assignMaterialProperties(*((*curIt).second), 1.);
403 if (confinedVolumes) {
405 std::span<Trk::TrackingVolume * const> volumes = confinedVolumes->
arrayObjects();
406 ATH_MSG_INFO(
"--> found : "<< volumes.size() <<
"confined TrackingVolumes");
408 for (
const auto & volume : volumes) {
411 assignLayerMaterialProperties(*volume, propSet);
420 ATH_MSG_INFO(
"========================================================================================= ");
423 #ifdef TRKDETDESCR_MEMUSAGE
424 m_memoryLogger.refresh(getpid());
425 ATH_MSG_INFO(
"[ memory usage ] Start building of material maps: " );
430 std::map< std::string, Trk::LayerMaterialMap* > layerMaterialMaps;
431 for (
auto& lmcIter : m_layerMaterialCreators ){
432 ATH_MSG_INFO(
"-> Creating material map '"<< lmcIter->layerMaterialName() <<
"' from creator "<< lmcIter.typeAndName() );
436 ATH_MSG_INFO( m_layerRecords.size() <<
" LayerMaterialRecords to be finalized for this material mapping run.");
439 for (
auto& lIter : m_layerRecords ) {
446 std::string
vName = eVolume ? (eVolume->
volumeName()) :
" BoundaryCollection ";
449 (lIter.second).finalizeRun(m_mapComposition);
451 if (!m_layerMaterialRecordAnalyser.empty() && m_layerMaterialRecordAnalyser->analyseLayerMaterial(*
layer, lIter.second).isFailure() )
454 bool analyse = (m_layerMaterialCreators.size() == m_layerMaterialAnalysers.size());
457 for (
auto& lmcIter : m_layerMaterialCreators ){
459 #ifdef TRKDETDESCR_MEMUSAGE
460 m_memoryLogger.refresh(getpid());
461 ATH_MSG_INFO(
"[ memory usage ] Before building the map for Layer "<< layerKey.
value() );
465 #ifdef TRKDETDESCR_MEMUSAGE
466 m_memoryLogger.refresh(getpid());
467 ATH_MSG_INFO(
"[ memory usage ] After building the map for Layer "<< layerKey.
value() );
473 (*layerMaterialMaps[lmcIter->layerMaterialName()])[layerKey.
value()] = lMaterial;
475 if (analyse && lMaterial && (m_layerMaterialAnalysers[ilmc]->analyseLayerMaterial(*
layer, *lMaterial)).isFailure() )
476 ATH_MSG_WARNING(
"Could not analyse created LayerMaterialProperties for layer "<< layerKey.
value() );
481 ATH_MSG_INFO(
"Finalize map synchronization and write the maps to the DetectorStore.");
483 for (
auto& ilmIter : layerMaterialMaps ){
485 if (m_mapComposition){
486 auto tElementTable = std::make_shared<Trk::ElementTable>(*m_elementTable);
487 ilmIter.second->updateElementTable(tElementTable);
488 if (ilmIter.second->elementTable()){
489 ATH_MSG_INFO(
"ElementTable for LayerMaterialMap '" << ilmIter.first <<
"' found and syncrhonized." );
494 if ( (
detStore()->record(ilmIter.second, ilmIter.first,
false)).isFailure()){
495 ATH_MSG_ERROR(
"Writing of LayerMaterialMap with name '" << ilmIter.first <<
"' was not successful." );
496 delete ilmIter.second;
497 }
else ATH_MSG_INFO(
"LayerMaterialMap: " << ilmIter.first <<
" written to the DetectorStore!" );
499 delete m_elementTable;
501 #ifdef TRKDETDESCR_MEMUSAGE
502 m_memoryLogger.refresh(getpid());
503 ATH_MSG_INFO(
"[ memory usage ] At the end of the material map creation.");
507 ATH_MSG_INFO(
"========================================================================================= " );
509 double unmapped = (m_unmapped+m_mapped) ?
double(m_unmapped)/
double(m_unmapped+m_mapped) : 0.;
511 ATH_MSG_INFO(
" -> Skipped (outisde) : " << m_skippedOutside );
512 ATH_MSG_INFO(
"========================================================================================= " );
514 return StatusCode::SUCCESS;
524 m_mappingVolume = trackingGeometry().trackingVolume(m_mappingVolumeName);
527 registerVolume(*trackingVolume, 0);
529 ATH_MSG_INFO(
"Add "<< m_layerRecords.size() <<
" confined volume layers to mapping setup.");
530 ATH_MSG_INFO(
"Add "<< trackingGeometry().numBoundaryLayers() <<
" boundary layers to mapping setup.");
533 for (
const auto bLayerIter : trackingGeometry().boundaryLayers())
534 insertLayerMaterialRecord(*(bLayerIter.first));
536 ATH_MSG_INFO(
"Map for "<< m_layerRecords.size() <<
" layers booked & prepared for mapping procedure");
538 return StatusCode::SUCCESS;
545 int sublevel = lvl+1;
549 std::cout <<
"TrackingVolume name: "<< tvol.
volumeName() << std::endl;
552 std::vector<const Trk::Layer*> volumeLayers;
556 if (confinedLayers) {
561 std::cout <<
"- found : "<<
layers.size() <<
"confined Layers"<< std::endl;
563 auto clIter =
layers.begin();
564 auto clIterE =
layers.end();
565 for ( ; clIter != clIterE; ++clIter ) {
567 const Amg::Vector3D& sReferencePoint = (*clIter)->surfaceRepresentation().globalReferencePoint();
568 bool insideMappingVolume = m_mappingVolume ? m_mappingVolume->inside(sReferencePoint) :
true;
569 if ((*clIter)->layerMaterialProperties() && insideMappingVolume)
570 volumeLayers.push_back((*clIter));
575 for (
auto& lIter : volumeLayers )
576 insertLayerMaterialRecord(*lIter);
580 if (confinedVolumes) {
581 std::span<Trk::TrackingVolume const * const> volumes = confinedVolumes->
arrayObjects();
585 std::cout <<
"- found : "<< volumes.size() <<
"confined TrackingVolumes"<< std::endl;
587 auto volumesIter = volumes.begin();
588 for (; volumesIter != volumes.end(); ++volumesIter)
590 registerVolume(**volumesIter, sublevel);
607 if (layerMaterialBinUtility){
610 layerMaterialBinUtility,
613 m_layerRecords[&lay] = lmr;
618 std::stringstream
msg;
619 msg <<
"Failed to get conditions data " << m_trackingGeometryReadKey.key() <<
".";
620 throw std::runtime_error(
msg.str());