ATLAS Offline Software
MaterialMapping.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // MaterialMapping.cxx, (c) ATLAS Detector software
8 
9 // Gaudi Units
10 #include "GaudiKernel/SystemOfUnits.h"
11 //TrkDetDescr Algs, Interfaces, Utils
12 #include "MaterialMapping.h"
16 // TrkGeometry
28 // TrkEvent
30 // Amg
32 
33 #ifdef TRKDETDESCR_MEMUSAGE
34 #include <unistd.h>
35 #endif
36 
37 Trk::MaterialMapping::MaterialMapping(const std::string& name, ISvcLocator* pSvcLocator)
38 : AthAlgorithm(name,pSvcLocator),
39  m_checkForEmptyHits(true),
40  m_mappingVolumeName("Atlas"),
41  m_mappingVolume(nullptr),
42  m_inputMaterialStepCollection("MaterialStepRecords"),
43  m_etaCutOff(6.0),
44  m_etaSide(0),
45  m_useLayerThickness(false),
46  m_associationType(1),
47  m_mapMaterial(true),
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.),
54  m_mapped(0),
55  m_unmapped(0),
56  m_skippedOutside(0),
57  m_layerMaterialScreenOutput(0)
58 #ifdef TRKDETDESCR_MEMUSAGE
59  ,m_memoryLogger()
60 #endif
61 {
62  // the name of the volume to map
63  declareProperty("MappingVolumeName" , m_mappingVolumeName);
64  // the extrapolation engine
65  declareProperty("CheckForEmptyHits" , m_checkForEmptyHits);
66  // general steering
67  declareProperty("EtaCutOff" , m_etaCutOff);
68  declareProperty("EtaSide" , m_etaSide);
69  // the toolhandle of the MaterialMapper to be used
70  declareProperty("MapMaterial" , m_mapMaterial);
71  // Composition related parameters
72  declareProperty("MapComposition" , m_mapComposition);
73  declareProperty("MinCompositionFraction" , m_minCompositionFraction);
74  // Steer the layer thickness
75  declareProperty("UseActualLayerThicknesss" , m_useLayerThickness);
76  // some job setup
77  declareProperty("MaterialAssociationType" , m_associationType);
78  declareProperty("InputMaterialStepCollection" , m_inputMaterialStepCollection);
79  declareProperty("InputElementTable" , m_inputEventElementTable);
80  // Output screen stuff
81  declareProperty("MaterialScreenOutputLevel" , m_layerMaterialScreenOutput);
82 
83 }
84 
86 = default;
87 
88 
90 {
91 
92  ATH_MSG_INFO("initialize()");
93 
94  ATH_CHECK( m_trackingGeometryReadKey.initialize(!m_trackingGeometryReadKey.key().empty()) );
95 
96  ATH_CHECK(m_extrapolationEngine.retrieve());
97 
98  if ( !m_materialMapper.empty() )
99  ATH_CHECK( m_materialMapper.retrieve() );
100 
101  if ( !m_layerMaterialRecordAnalyser.empty() )
102  ATH_CHECK( m_layerMaterialRecordAnalyser.retrieve() );
103 
104  if ( !m_layerMaterialCreators.empty() )
105  ATH_CHECK( m_layerMaterialCreators.retrieve() );
106 
107  if ( !m_layerMaterialAnalysers.empty() )
108  ATH_CHECK( m_layerMaterialAnalysers.retrieve() );
109 
110  ATH_CHECK( m_inputMaterialStepCollection.initialize() );
111  ATH_CHECK( m_inputEventElementTable.initialize() );
112 
113  return StatusCode::SUCCESS;
114 }
115 
116 
118 {
119  ATH_MSG_VERBOSE("MaterialMapping execute() start");
120 
121  // ------------------------------- get the trackingGeometry at first place
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.");
126  return retrieveCode;
127  }
128  }
129  if (m_mappingVolume)
130  ATH_MSG_VERBOSE("Mapping volume correctly retrieved from tracking geometry");
131 
132 
133  SG::ReadHandle<MaterialStepCollection> materialStepCollection(m_inputMaterialStepCollection);
134 
135  // --------- prepare the element table ---------------------------------------------------
136 
137  if (m_mapComposition) {
138  SG::ReadHandle<Trk::ElementTable> eTableEvent(m_inputEventElementTable);
139  (*m_elementTable) += (*eTableEvent); // accummulate the table
140  }
141 
142 
143  // event parameters - associated asteps, and layers hit per event
144  int associatedSteps = 0;
145  m_accumulatedMaterialXX0 = 0.;
146  m_accumulatedRhoS = 0.;
147  m_layersRecordedPerEvent.clear();
148  // clearing the recorded layers per event
149  if (materialStepCollection.isValid() && !materialStepCollection->empty()){
150 
151  // get the number of material steps
152  size_t materialSteps = materialStepCollection->size();
153  ATH_MSG_DEBUG("[+] Successfully read "<< materialSteps << " geantino steps");
154 
155  // create a direction out of the last material step
156  double dirx = (*materialStepCollection)[materialSteps-1]->hitX();
157  double diry = (*materialStepCollection)[materialSteps-1]->hitY();
158  double dirz = (*materialStepCollection)[materialSteps-1]->hitZ();
159  Amg::Vector3D direction = Amg::Vector3D(dirx,diry,dirz).unit();
160 
161  double eta = direction.eta();
162  // skip the event if the eta cut is not met
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;
166  }
167 
168  // now propagate through the full detector and collect the layers
169  Trk::NeutralCurvilinearParameters ncP(Amg::Vector3D(0.,0.,0.), direction, 0.);
170  // create a neutral extrapolation cell
172  ecc.navigationCurvilinear = false;
176 
177  // let's extrapolate through the detector and remember which layers (with material) should have been hit
178  std::vector< std::pair<const Trk::Layer*, Amg::Vector3D> > layersAndHits;
179  // call the extrapolation engine
180  Trk::ExtrapolationCode eCode = m_extrapolationEngine->extrapolate(ecc);
181  // end the parameters if there
182  if (eCode.isSuccess()){
183  // name of passive surfaces found
184  size_t nLayersHit = ecc.extrapolationSteps.size();
185  ATH_MSG_VERBOSE("[+] Extrapolation to layers did succeed and found " << nLayersHit << " layers.");
186  // reserve the size of the vectors
187  layersAndHits.reserve(nLayersHit);
188  // for screen output
189  size_t ilayer = 0;
190  // find all the intersected material - remember the last parameters
191  const Trk::NeutralParameters* parameters = nullptr;
192  // loop over the collected information
193  for (auto& es : ecc.extrapolationSteps){
194  // continue if we have parameters
195  parameters = es.parameters;
196  if (parameters){
197  const Trk::Surface& pSurface = parameters->associatedSurface();
198  // get the surface with associated layer (that has material)
199  ATH_MSG_VERBOSE("[L] Testing layer with associatedLayer() " << pSurface.associatedLayer() << " and materialLayer() " << pSurface.materialLayer() );
200 
201  if ((pSurface.associatedLayer() &&
202  pSurface.associatedLayer()->layerMaterialProperties()) ||
203  pSurface.materialLayer()) {
204  // material layer
205 
206  const Trk::Layer* mLayer = pSurface.materialLayer()
207  ? pSurface.materialLayer()
208  : pSurface.associatedLayer();
209  // record that one
210  std::pair<const Trk::Layer*, Amg::Vector3D> layerHitPair(
211  mLayer, parameters->position());
212  ATH_MSG_VERBOSE("[L] Layer "
213  << ++ilayer << " with index "
214  << mLayer->layerIndex().value()
215  << " hit at "
216  << Amg::toString(parameters->position()));
217  layersAndHits.push_back(layerHitPair);
218  }
219  delete parameters;
220  }
221  }
222  // cleanup of the final hits
223  if (ecc.endParameters != parameters) delete ecc.endParameters;
224 
225  // we have no layers and Hits
226  if (layersAndHits.empty()){
227  ATH_MSG_VERBOSE("[!] No Layer was intersected - skipping.");
228  return StatusCode::SUCCESS;
229  }
230 
231  // layers are ordered, hence you can move the starting point along
232  size_t currentLayer = 0;
233  // loop through hits and find the closest layer, the start point moves outwards as we go
234  for ( const Trk::MaterialStep* step : *materialStepCollection ) {
235  // verbose output
236  ATH_MSG_VERBOSE("[L] starting from layer " << currentLayer << " from layer collection for this step.");
237  // step length and position
238  double t = step->steplength();
239  Amg::Vector3D pos(step->hitX(), step->hitY(), step->hitZ());
240  // skip if :
241  // -- 0) no mapping volume exists
242  // -- 1) outside the mapping volume
243  // -- 2) outside the eta acceptance
244  if (!m_mappingVolume || !(m_mappingVolume->inside(pos)) || fabs(pos.eta()) > m_etaCutOff ){
245  ++m_skippedOutside;
246  continue;
247  }
248  // now find the closest layer
249  // (a) if the currentLayer is the last layer and the hit is still inside -> assign
250  if (currentLayer < nLayersHit-1) {
251  // search through the layers - this is the reference distance for projection
252  double currentDistance = (pos-layersAndHits[currentLayer].second).mag();
253  ATH_MSG_VERBOSE("- current distance is " << currentDistance << " from " << Amg::toString(pos) << " and " << Amg::toString(layersAndHits[currentLayer].second) );
254  for (size_t testLayer = (currentLayer+1); testLayer < nLayersHit; ++testLayer){
255  // calculate teh distance to the testLayer
256  double testDistance = (pos-layersAndHits[testLayer].second).mag();
257  ATH_MSG_VERBOSE("[L] Testing layer " << testLayer << " from layer collection for this step.");
258  ATH_MSG_VERBOSE("- test distance is " << testDistance << " from " << Amg::toString(pos) << " and " << Amg::toString(layersAndHits[testLayer].second) );
259  if ( testDistance < currentDistance ){
260  // screen output
261  ATH_MSG_VERBOSE("[L] Skipping over to current layer " << testLayer << " because " << testDistance << " < " << currentDistance);
262  // the test distance did shrink - update currentLayer
263  currentLayer = testLayer;
264  currentDistance = testDistance;
265  } else {
266  // stick to the layer you have
267  break;
268  }
269  }
270  }
271  // the currentLayer *should* be correct now
272  const Trk::Layer* assignedLayer = layersAndHits[currentLayer].first;
273  Amg::Vector3D assignedPosition = layersAndHits[currentLayer].second;
274  // associate the hit
275  // (1) count it
276  ++associatedSteps;
277  // (2) associate it
278  associateHit(*assignedLayer, pos, assignedPosition, t, step->fullMaterial());
279  } // loop over material Steps
280 
281  // check for the empty hits - they need to be taken into account
282  ATH_MSG_VERBOSE("Found " << layersAndHits.size() << " intersected layers - while having " << m_layersRecordedPerEvent.size() << " recorded ones.");
283 
284  // now - cross-chek if you have additional layers
285  for ( auto& lhp : layersAndHits){
286  // check if you find the layer int he already done record-map : not found - we need to do an empty hit scaling
287  if (m_layersRecordedPerEvent.find(lhp.first) == m_layersRecordedPerEvent.end()){
288  // try to find the layer material record
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.");
293  } else
294  ATH_MSG_WARNING("- no Layer found in the associated map! Should not happen.");
295  }
296  }
297 
298  // check whether the event was good for at least one hit
299  if (associatedSteps) {
300  ATH_MSG_VERBOSE("There are associated steps, need to call finalizeEvent() & record to the MaterialMapper.");
301  // finalize the event --------------------- Layers ---------------------------------------------
302  for (auto& lRecord : m_layerRecords ) {
303  // associated material
304  Trk::AssociatedMaterial* assMatHit = lRecord.second.finalizeEvent((*lRecord.first));
305  // record the full layer hit
306  if (assMatHit && !m_materialMapper.empty()) m_materialMapper->recordLayerHit(*assMatHit, true);
307  delete assMatHit;
308  // call the material mapper finalize method
309  ATH_MSG_VERBOSE("Calling finalizeEvent on the MaterialMapper ...");
310  }
311  } // the event had at least one associated hit
312 
313  } // end of eCode.success : needed for new mapping schema
314 
315  } // material steps existed
316 
317 
318  return StatusCode::SUCCESS;
319 }
320 
321 
322 bool Trk::MaterialMapping::associateHit( const Trk::Layer& associatedLayer,
323  const Amg::Vector3D& pos,
324  const Amg::Vector3D& positionOnLayer,
325  double stepl,
326  const Trk::Material& mat)
327 {
328  // get the intersection with the layer
329  ++m_mapped;
330  // get the layer
331  const Trk::Layer* layer = &associatedLayer;
332 
333  // get the associated volume
334  const Trk::TrackingVolume* associatedVolume = trackingGeometry().lowestTrackingVolume(pos);
335 
336  // try to find the layer material record
337  auto clIter = m_layerRecords.find(layer);
338  if (clIter != m_layerRecords.end() ){
339  // remember that you actually hit this layer
340  m_layersRecordedPerEvent[layer] = true;
341  // LayerMaterialRecord found, add the hit
342  (*clIter).second.associateGeantinoHit(positionOnLayer, stepl, mat);
343  ATH_MSG_VERBOSE("- associate Geantino Information at intersection [r/z] = " << positionOnLayer.perp() << "/"<< positionOnLayer.z() );
344  ATH_MSG_VERBOSE(" mapping distance = " << (pos-positionOnLayer).mag() );
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);
351  ATH_MSG_VERBOSE("- to layer with index "<< layer->layerIndex().value() << " from '"<< associatedVolume->volumeName() << "'.");
352  // record the plain information w/o correction & intersection
353  if (m_mapMaterial && !m_materialMapper.empty()){
354  Trk::AssociatedMaterial am(pos, stepl, mat, 1., associatedVolume, layer);
355  m_materialMapper->recordMaterialHit(am, positionOnLayer);
356  ATH_MSG_VERBOSE(" - associated material produced as : " << am);
357  }
358 
359  } else if (layer) {
360  ATH_MSG_WARNING("- associate hit - the layer with index " << layer->layerIndex().value() << " was not found - should not happen!");
361  }
362  // return
363  return true;
364 }
365 
366 
368  Trk::LayerMaterialMap* propSet)
369 {
370 
371  if (!propSet) return;
372 
373  ATH_MSG_INFO("Processing TrackingVolume: "<< tvol.volumeName() );
374 
375  // ----------------------------------- loop over confined layers ------------------------------------------
376  Trk::BinnedArray< Trk::Layer >* confinedLayers = tvol.confinedLayers();
377  if (confinedLayers) {
378  // get the objects in a vector-like format
379  std::span<Trk::Layer * const> layers = confinedLayers->arrayObjects();
380  ATH_MSG_INFO("--> found : "<< layers.size() << "confined Layers");
381  // the iterator over the vector
382  // loop over layers
383  for (Trk::Layer* layer : layers) {
384  // assign the material and output
385  if (layer && (*layer).layerIndex().value() ) {
386  ATH_MSG_INFO(" > LayerIndex: "<< (*layer).layerIndex() );
387  // set the material!
388  if (propSet) {
389  // find the pair
390  auto curIt = propSet->find((*layer).layerIndex());
391  if (curIt != propSet->end()) {
392  ATH_MSG_INFO("LayerMaterial assigned for Layer with index: "<< (*layer).layerIndex() );
393  // set it to the layer
394  layer->assignMaterialProperties(*((*curIt).second), 1.);
395  }
396  }
397  }
398  }
399  }
400 
401  // ----------------------------------- loop over confined volumes -----------------------------
403  if (confinedVolumes) {
404  // get the objects in a vector-like format
405  std::span<Trk::TrackingVolume * const> volumes = confinedVolumes->arrayObjects();
406  ATH_MSG_INFO("--> found : "<< volumes.size() << "confined TrackingVolumes");
407  // loop over volumes
408  for (const auto & volume : volumes) {
409  // assing the material and output
410  if (volume)
411  assignLayerMaterialProperties(*volume, propSet); // call itself recursively
412  }
413  }
414 }
415 
416 
418 {
419 
420  ATH_MSG_INFO("========================================================================================= ");
421  ATH_MSG_INFO("finalize() starts ...");
422 
423 #ifdef TRKDETDESCR_MEMUSAGE
424  m_memoryLogger.refresh(getpid());
425  ATH_MSG_INFO("[ memory usage ] Start building of material maps: " );
426  ATH_MSG_INFO( m_memoryLogger );
427 #endif
428 
429  // create a dedicated LayerMaterialMap by layerMaterialCreator;
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() );
433  layerMaterialMaps[lmcIter->layerMaterialName()] = new Trk::LayerMaterialMap();
434  }
435 
436  ATH_MSG_INFO( m_layerRecords.size() << " LayerMaterialRecords to be finalized for this material mapping run.");
437 
438  // loop over the layers and output the stuff --- fill the associatedLayerMaterialProperties
439  for ( auto& lIter : m_layerRecords ) {
440  // Get the key map, the layer & the volume name
441  const Trk::Layer* layer = lIter.first;
442  Trk::LayerIndex layerKey = layer->layerIndex();
443  // get the enclosing tracking volume
444  const Trk::TrackingVolume* eVolume = layer->enclosingTrackingVolume();
445  // assign the string
446  std::string vName = eVolume ? (eVolume->volumeName()) : " BoundaryCollection ";
447  ATH_MSG_INFO("Finalize MaterialAssociation for Layer "<< layerKey.value() << " in " << vName );
448  // normalize - use m_finalizeRunDebug
449  (lIter.second).finalizeRun(m_mapComposition);
450  // output the material to the analyser if registered
451  if (!m_layerMaterialRecordAnalyser.empty() && m_layerMaterialRecordAnalyser->analyseLayerMaterial(*layer, lIter.second).isFailure() )
452  ATH_MSG_WARNING("Could not analyse the LayerMaterialRecord for layer "<< layerKey.value() );
453  // check if we have analysers per creator
454  bool analyse = (m_layerMaterialCreators.size() == m_layerMaterialAnalysers.size());
455  // and now use the creators to make the maps out of the LayerMaterialRecord
456  size_t ilmc = 0;
457  for ( auto& lmcIter : m_layerMaterialCreators ){
458  // call the creator and register in the according map
459 #ifdef TRKDETDESCR_MEMUSAGE
460  m_memoryLogger.refresh(getpid());
461  ATH_MSG_INFO("[ memory usage ] Before building the map for Layer "<< layerKey.value() );
462  ATH_MSG_INFO( m_memoryLogger );
463 #endif
464  const Trk::LayerMaterialProperties* lMaterial = lmcIter->createLayerMaterial(lIter.second);
465 #ifdef TRKDETDESCR_MEMUSAGE
466  m_memoryLogger.refresh(getpid());
467  ATH_MSG_INFO("[ memory usage ] After building the map for Layer "<< layerKey.value() );
468  ATH_MSG_INFO( m_memoryLogger );
469 #endif
470  if (lMaterial)
471  ATH_MSG_VERBOSE("LayerMaterial map created as "<< *lMaterial );
472  // insert the created map for the given layer
473  (*layerMaterialMaps[lmcIter->layerMaterialName()])[layerKey.value()] = lMaterial;
474  // analyse the it if configured
475  if (analyse && lMaterial && (m_layerMaterialAnalysers[ilmc]->analyseLayerMaterial(*layer, *lMaterial)).isFailure() )
476  ATH_MSG_WARNING("Could not analyse created LayerMaterialProperties for layer "<< layerKey.value() );
477  ++ilmc;
478  }
479  }
480 
481  ATH_MSG_INFO("Finalize map synchronization and write the maps to the DetectorStore.");
482 
483  for (auto& ilmIter : layerMaterialMaps ){
484  // elementTable handling - if existent
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." );
490  ATH_MSG_INFO( *(ilmIter.second->elementTable()) );
491  }
492  }
493  // detector store writing
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!" );
498  }
499  delete m_elementTable;
500 
501 #ifdef TRKDETDESCR_MEMUSAGE
502  m_memoryLogger.refresh(getpid());
503  ATH_MSG_INFO( "[ memory usage ] At the end of the material map creation.");
504  ATH_MSG_INFO( m_memoryLogger );
505 #endif
506 
507  ATH_MSG_INFO( "========================================================================================= " );
508  ATH_MSG_INFO( " -> Total mapped hits : " << m_mapped );
509  double unmapped = (m_unmapped+m_mapped) ? double(m_unmapped)/double(m_unmapped+m_mapped) : 0.;
510  ATH_MSG_INFO( " -> Total (rel.) unmapped hits : " << m_unmapped << " (" << unmapped << ")" );
511  ATH_MSG_INFO( " -> Skipped (outisde) : " << m_skippedOutside );
512  ATH_MSG_INFO( "========================================================================================= " );
513  ATH_MSG_INFO( "finalize() successful");
514  return StatusCode::SUCCESS;
515 }
516 
517 
519 {
520  // either get a string volume or the highest one
521  const Trk::TrackingVolume* trackingVolume = trackingGeometry().highestTrackingVolume();
522 
523  // prepare the mapping volume
524  m_mappingVolume = trackingGeometry().trackingVolume(m_mappingVolumeName);
525 
526  // register the confined layers from the TrackingVolume
527  registerVolume(*trackingVolume, 0);
528 
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.");
531 
532  // register the layers from boundary surfaces
533  for (const auto bLayerIter : trackingGeometry().boundaryLayers())
534  insertLayerMaterialRecord(*(bLayerIter.first));
535 
536  ATH_MSG_INFO("Map for "<< m_layerRecords.size() << " layers booked & prepared for mapping procedure");
537 
538  return StatusCode::SUCCESS;
539 
540 }
541 
542 
544 {
545  int sublevel = lvl+1;
546 
547  for (int indent=0; indent<sublevel; ++indent)
548  std::cout << " ";
549  std::cout << "TrackingVolume name: "<< tvol.volumeName() << std::endl;
550 
551  // all those to be processed
552  std::vector<const Trk::Layer*> volumeLayers;
553 
554  // collect all material layers that have layerMaterial
555  const Trk::BinnedArray< Trk::Layer >* confinedLayers = tvol.confinedLayers();
556  if (confinedLayers) {
557  // this go ahead with the layers
558  std::span<Trk::Layer const * const> layers = confinedLayers->arrayObjects();
559  for (int indent=0; indent<sublevel; ++indent)
560  std::cout << " ";
561  std::cout << "- found : "<< layers.size() << "confined Layers"<< std::endl;
562  // loop over and fill them
563  auto clIter = layers.begin();
564  auto clIterE = layers.end();
565  for ( ; clIter != clIterE; ++clIter ) {
566  // only take layers with MaterialProperties defined and which are within the mapping volume
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));
571  }
572  }
573 
574  // now create LayerMaterialRecords for all
575  for ( auto& lIter : volumeLayers )
576  insertLayerMaterialRecord(*lIter);
577 
578  // step dopwn the navigation tree to reach the confined volumes
579  const Trk::BinnedArray<Trk::TrackingVolume >* confinedVolumes = tvol.confinedVolumes();
580  if (confinedVolumes) {
581  std::span<Trk::TrackingVolume const * const> volumes = confinedVolumes->arrayObjects();
582 
583  for (int indent=0; indent<sublevel; ++indent)
584  std::cout << " ";
585  std::cout << "- found : "<< volumes.size() << "confined TrackingVolumes"<< std::endl;
586  // loop over the confined volumes
587  auto volumesIter = volumes.begin();
588  for (; volumesIter != volumes.end(); ++volumesIter)
589  if (*volumesIter) {
590  registerVolume(**volumesIter, sublevel);
591  }
592  }
593 
594 }
595 
597  // first occurrance, create a new LayerMaterialRecord
598  // - get the bin utility for the binned material (if necessary)
599  // - get the material first
600  const Trk::LayerMaterialProperties* layerMaterialProperties = lay.layerMaterialProperties();
601  // - dynamic cast to the BinnedLayerMaterial
602  const Trk::BinnedLayerMaterial* layerBinnedMaterial
603  = dynamic_cast<const Trk::BinnedLayerMaterial*>(layerMaterialProperties);
604  // get the binned array
605  const Trk::BinUtility* layerMaterialBinUtility = (layerBinnedMaterial) ? layerBinnedMaterial->binUtility() : nullptr;
606  // now fill the layer material record
607  if (layerMaterialBinUtility){
608  // create a new Layer Material record in the map
609  Trk::LayerMaterialRecord lmr((m_useLayerThickness ? lay.thickness() : 1.),
610  layerMaterialBinUtility,
611  (Trk::MaterialAssociationType)m_associationType);
612  // and fill it into the map
613  m_layerRecords[&lay] = lmr;
614  }
615 }
616 
618  std::stringstream msg;
619  msg << "Failed to get conditions data " << m_trackingGeometryReadKey.key() << ".";
620  throw std::runtime_error(msg.str());
621 }
622 
623 
Trk::ExtrapolationMode::CollectPassive
@ CollectPassive
Definition: ExtrapolationCell.h:57
MaterialMapping.h
Trk::Surface::materialLayer
const Trk::MaterialLayer * materialLayer() const
return the material Layer
Trk::ExtrapolationCell::extrapolationSteps
std::vector< ExtrapolationStep< T > > extrapolationSteps
parameters on sensitive detector elements
Definition: ExtrapolationCell.h:292
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
Trk::ExtrapolationCell::addConfigurationMode
void addConfigurationMode(ExtrapolationMode::eMode em)
add a configuration mode
Definition: ExtrapolationCell.h:345
MaterialProperties.h
Trk::ExtrapolationMode::CollectBoundary
@ CollectBoundary
Definition: ExtrapolationCell.h:58
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:67
Trk::MaterialMapping::m_inputEventElementTable
SG::ReadHandleKey< Trk::ElementTable > m_inputEventElementTable
input event table
Definition: MaterialMapping.h:167
mat
GeoMaterial * mat
Definition: LArDetectorConstructionTBEC.cxx:55
Trk::MaterialMapping::registerVolume
void registerVolume(const Trk::TrackingVolume &tvol, int lvl)
Output information with Level.
Definition: MaterialMapping.cxx:543
BinUtility.h
module_driven_slicing.layers
layers
Definition: module_driven_slicing.py:113
python.SystemOfUnits.second
float second
Definition: SystemOfUnits.py:135
Trk::MaterialMapping::m_etaCutOff
double m_etaCutOff
general steering
Definition: MaterialMapping.h:131
HomogeneousLayerMaterial.h
Trk::MaterialMapping::handleTrackingGeometry
StatusCode handleTrackingGeometry()
Retrieve the TrackingGeometry and its informations.
Definition: MaterialMapping.cxx:518
NeutralParameters.h
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
python.compareNtuple.vName
vName
Definition: compareNtuple.py:23
Trk::MaterialMapping::throwFailedToGetTrackingGeometry
void throwFailedToGetTrackingGeometry() const
Definition: MaterialMapping.cxx:617
Trk::LayerMaterialMap
Definition: LayerMaterialMap.h:32
Trk::MaterialMapping::m_minCompositionFraction
double m_minCompositionFraction
minimal fraction to be accounted for the composition recording
Definition: MaterialMapping.h:162
Trk::TrackingVolume::confinedLayers
const LayerArray * confinedLayers() const
Return the subLayer array.
Trk::MaterialStep
Definition: MaterialStep.h:34
GeometryStatics.h
Trk::MaterialMapping::m_associationType
int m_associationType
Definition: MaterialMapping.h:134
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition: GeoPrimitivesToStringConverter.h:40
Trk::MaterialMapping::assignLayerMaterialProperties
void assignLayerMaterialProperties(Trk::TrackingVolume &tvol, Trk::LayerMaterialMap *propSet)
create the LayerMaterialRecord *‍/
Definition: MaterialMapping.cxx:367
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
geometry_dat_to_json.indent
indent
Definition: geometry_dat_to_json.py:37
Trk::MaterialMapping::m_etaSide
int m_etaSide
needed for debugging: -1 negative | 0 all | 1 positive
Definition: MaterialMapping.h:132
Trk::LayerIndex
Definition: LayerIndex.h:37
LayerMaterialRecord.h
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
Trk::MaterialMapping::m_mapComposition
bool m_mapComposition
map the composition of the material
Definition: MaterialMapping.h:161
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
Definition: AthCommonDataStore.h:145
TRT::Hit::layer
@ layer
Definition: HitInfo.h:79
Trk::ExtrapolationCode::isSuccess
bool isSuccess() const
return success
Definition: ExtrapolationCell.h:153
Trk::LayerMaterialProperties
Definition: LayerMaterialProperties.h:62
Trk::ExtrapolationCode
Definition: ExtrapolationCell.h:105
Trk::ExtrapolationCell::endParameters
T * endParameters
by pointer - are newly created and can be optionally 0
Definition: ExtrapolationCell.h:238
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
BinnedLayerMaterial.h
Trk::ParametersBase
Definition: ParametersBase.h:55
Trk::CurvilinearParametersT
Definition: CurvilinearParametersT.h:48
Trk::BinnedLayerMaterial
Definition: BinnedLayerMaterial.h:33
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
Trk::MaterialMapping::m_mappingVolumeName
std::string m_mappingVolumeName
Definition: MaterialMapping.h:123
Trk::LayerIndex::value
int value() const
layerIndex expressed in an integer
Definition: LayerIndex.h:71
LayerIndex.h
Trk::BinnedArray::arrayObjects
virtual std::span< T *const > arrayObjects()=0
Return all objects of the Array non-const we can still modify the T.
AthAlgorithm
Definition: AthAlgorithm.h:47
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
Trk::MaterialMapping::~MaterialMapping
~MaterialMapping()
Default Destructor.
Trk::AssociatedMaterial
Definition: AssociatedMaterial.h:33
Trk::BinUtility
Definition: BinUtility.h:39
Trk::MaterialMapping::initialize
StatusCode initialize()
standard Athena-Algorithm method
Definition: MaterialMapping.cxx:89
Trk::LayerMaterialRecord
Definition: LayerMaterialRecord.h:42
python.PyKernel.detStore
detStore
Definition: PyKernel.py:41
Trk::TrackingVolume::volumeName
const std::string & volumeName() const
Returns the VolumeName - for debug reason, might be depreciated later.
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
MaterialLayer.h
Trk::MaterialMapping::execute
StatusCode execute()
standard Athena-Algorithm method
Definition: MaterialMapping.cxx:117
Trk::Layer::thickness
double thickness() const
Return the Thickness of the Layer.
LayerMaterialMap.h
Trk::MaterialAssociationType
MaterialAssociationType
Definition: MaterialAssociationType.h:13
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:16
Trk::MaterialMapping::m_useLayerThickness
bool m_useLayerThickness
use the actual layer thickness
Definition: MaterialMapping.h:133
Trk::MaterialMapping::m_mapMaterial
bool m_mapMaterial
Mapper and Inspector.
Definition: MaterialMapping.h:156
Trk::Surface::associatedLayer
const Trk::Layer * associatedLayer() const
return the associated Layer
Trk::MaterialMapping::associateHit
bool associateHit(const Trk::Layer &tvol, const Amg::Vector3D &pos, const Amg::Vector3D &layerHitPosition, double stepl, const Trk::Material &mat)
Associate the Step to the Layer.
Definition: MaterialMapping.cxx:322
Trk::BinnedLayerMaterial::binUtility
virtual const BinUtility * binUtility() const override
Return the BinUtility.
Definition: BinnedLayerMaterial.h:121
TrackingVolume.h
Trk::ExtrapolationCell
Definition: ExtrapolationCell.h:231
Trk::MaterialMapping::m_checkForEmptyHits
bool m_checkForEmptyHits
use extrapoaltion engine to check for empty hits
Definition: MaterialMapping.h:115
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Trk::TrackingVolume::confinedVolumes
const TrackingVolumeArray * confinedVolumes() const
Return the subLayer array.
Trk::MaterialMapping::MaterialMapping
MaterialMapping(const std::string &name, ISvcLocator *pSvcLocator)
Standard Athena-Algorithm Constructor.
Definition: MaterialMapping.cxx:37
Trk::Layer::layerMaterialProperties
const LayerMaterialProperties * layerMaterialProperties() const
getting the LayerMaterialProperties including full/pre/post update
Trk::ExtrapolationCell::navigationCurvilinear
bool navigationCurvilinear
stay in curvilinear parameters where possible, default : true
Definition: ExtrapolationCell.h:284
GeoPrimitivesToStringConverter.h
LArCellBinning.step
step
Definition: LArCellBinning.py:158
Trk::MaterialMapping::insertLayerMaterialRecord
void insertLayerMaterialRecord(const Trk::Layer &lay)
Definition: MaterialMapping.cxx:596
AssociatedMaterial.h
physics_parameters.parameters
parameters
Definition: physics_parameters.py:144
Trk::Material
Definition: Material.h:117
Trk::MaterialMapping::m_inputMaterialStepCollection
SG::ReadHandleKey< MaterialStepCollection > m_inputMaterialStepCollection
output / input steering
Definition: MaterialMapping.h:127
Trk::MaterialMapping::finalize
StatusCode finalize()
standard Athena-Algorithm method
Definition: MaterialMapping.cxx:417
Trk::BinnedArray
Definition: BinnedArray.h:36
Trk::Surface
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:79
Trk::ExtrapolationMode::StopAtBoundary
@ StopAtBoundary
Definition: ExtrapolationCell.h:55
Trk::TrackingVolume
Definition: TrackingVolume.h:119
python.AutoConfigFlags.msg
msg
Definition: AutoConfigFlags.py:7
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
LayerMaterialProperties.h
Trk::MaterialMapping::m_layerMaterialScreenOutput
int m_layerMaterialScreenOutput
Definition: MaterialMapping.h:183
HitType::unmapped
@ unmapped
Trk::Layer
Definition: Layer.h:72
Trk::Layer::layerIndex
const LayerIndex & layerIndex() const
get the layerIndex
MaterialStep.h
CompressedLayerMaterial.h