10 #include "GaudiKernel/SystemOfUnits.h"
32 #ifdef TRKDETDESCR_MEMUSAGE
38 m_checkForEmptyHits(true),
39 m_mappingVolumeName(
"Atlas"),
40 m_mappingVolume(nullptr),
41 m_inputMaterialStepCollection(
"MaterialStepRecords"),
44 m_useLayerThickness(false),
47 m_mapComposition(false),
48 m_minCompositionFraction(0.005),
49 m_elementTable(nullptr),
50 m_inputEventElementTable(
"ElementTable"),
51 m_accumulatedMaterialXX0(0.),
52 m_accumulatedRhoS(0.),
56 m_layerMaterialScreenOutput(0)
57 #ifdef TRKDETDESCR_MEMUSAGE
93 ATH_CHECK( m_trackingGeometryReadKey.initialize(!m_trackingGeometryReadKey.key().empty()) );
95 ATH_CHECK(m_extrapolationEngine.retrieve());
97 if ( !m_materialMapper.empty() )
100 if ( !m_layerMaterialRecordAnalyser.empty() )
101 ATH_CHECK( m_layerMaterialRecordAnalyser.retrieve() );
103 if ( !m_layerMaterialCreators.empty() )
104 ATH_CHECK( m_layerMaterialCreators.retrieve() );
106 if ( !m_layerMaterialAnalysers.empty() )
107 ATH_CHECK( m_layerMaterialAnalysers.retrieve() );
109 ATH_CHECK( m_inputMaterialStepCollection.initialize() );
110 ATH_CHECK( m_inputEventElementTable.initialize() );
112 return StatusCode::SUCCESS;
121 if (!m_mappingVolume) {
122 StatusCode retrieveCode = handleTrackingGeometry();
123 if (retrieveCode.isFailure()){
124 ATH_MSG_INFO(
"Could not retrieve mapping volume from tracking geometry. Exiting.");
129 ATH_MSG_VERBOSE(
"Mapping volume correctly retrieved from tracking geometry");
136 if (m_mapComposition) {
138 (*m_elementTable) += (*eTableEvent);
143 int associatedSteps = 0;
144 m_accumulatedMaterialXX0 = 0.;
145 m_accumulatedRhoS = 0.;
146 m_layersRecordedPerEvent.clear();
148 if (materialStepCollection.
isValid() && !materialStepCollection->empty()){
151 size_t materialSteps = materialStepCollection->size();
152 ATH_MSG_DEBUG(
"[+] Successfully read "<< materialSteps <<
" geantino steps");
155 double dirx = (*materialStepCollection)[materialSteps-1]->hitX();
156 double diry = (*materialStepCollection)[materialSteps-1]->hitY();
157 double dirz = (*materialStepCollection)[materialSteps-1]->hitZ();
160 double eta = direction.eta();
162 if ( fabs(
eta) > m_etaCutOff || (m_etaSide && m_etaSide*
eta < 0.) ) {
163 ATH_MSG_VERBOSE(
"[-] Event is outside eta acceptance of " << m_etaCutOff <<
". Skipping it.");
164 return StatusCode::SUCCESS;
177 std::vector< std::pair<const Trk::Layer*, Amg::Vector3D> > layersAndHits;
184 ATH_MSG_VERBOSE(
"[+] Extrapolation to layers did succeed and found " << nLayersHit <<
" layers.");
186 layersAndHits.reserve(nLayersHit);
204 std::pair<const Trk::Layer*, Amg::Vector3D> layerHitPair(mLayer,
parameters->position());
206 layersAndHits.push_back(layerHitPair);
215 if (layersAndHits.empty()){
217 return StatusCode::SUCCESS;
221 size_t currentLayer = 0;
225 ATH_MSG_VERBOSE(
"[L] starting from layer " << currentLayer <<
" from layer collection for this step.");
227 double t =
step->steplength();
233 if (!m_mappingVolume || !(m_mappingVolume->inside(
pos)) || fabs(
pos.eta()) > m_etaCutOff ){
239 if (currentLayer < nLayersHit-1) {
241 double currentDistance = (
pos-layersAndHits[currentLayer].second).
mag();
243 for (
size_t testLayer = (currentLayer+1); testLayer < nLayersHit; ++testLayer){
245 double testDistance = (
pos-layersAndHits[testLayer].second).
mag();
246 ATH_MSG_VERBOSE(
"[L] Testing layer " << testLayer <<
" from layer collection for this step.");
248 if ( testDistance < currentDistance ){
250 ATH_MSG_VERBOSE(
"[L] Skipping over to current layer " << testLayer <<
" because " << testDistance <<
" < " << currentDistance);
252 currentLayer = testLayer;
253 currentDistance = testDistance;
261 const Trk::Layer* assignedLayer = layersAndHits[currentLayer].first;
262 Amg::Vector3D assignedPosition = layersAndHits[currentLayer].second;
267 associateHit(*assignedLayer,
pos, assignedPosition,
t,
step->fullMaterial());
271 ATH_MSG_VERBOSE(
"Found " << layersAndHits.size() <<
" intersected layers - while having " << m_layersRecordedPerEvent.size() <<
" recorded ones.");
274 for (
auto& lhp : layersAndHits){
276 if (m_layersRecordedPerEvent.find(lhp.first) == m_layersRecordedPerEvent.end()){
278 auto clIter = m_layerRecords.find(lhp.first);
279 if (clIter != m_layerRecords.end() ){
280 (*clIter).second.associateEmptyHit(lhp.second);
281 ATH_MSG_VERBOSE(
"- to layer with index "<< lhp.first->layerIndex().value() <<
" with empty hit detected.");
283 ATH_MSG_WARNING(
"- no Layer found in the associated map! Should not happen.");
288 if (associatedSteps) {
289 ATH_MSG_VERBOSE(
"There are associated steps, need to call finalizeEvent() & record to the MaterialMapper.");
291 for (
auto& lRecord : m_layerRecords ) {
295 if (assMatHit && !m_materialMapper.empty()) m_materialMapper->recordLayerHit(*assMatHit,
true);
307 return StatusCode::SUCCESS;
326 auto clIter = m_layerRecords.find(
layer);
327 if (clIter != m_layerRecords.end() ){
329 m_layersRecordedPerEvent[
layer] =
true;
331 (*clIter).second.associateGeantinoHit(positionOnLayer, stepl,
mat);
332 ATH_MSG_VERBOSE(
"- associate Geantino Information at intersection [r/z] = " << positionOnLayer.perp() <<
"/"<< positionOnLayer.z() );
334 ATH_MSG_VERBOSE(
"- associate Geantino Information ( s, s/x0 , x0 , l0, a, z, rho ) = "
335 << stepl <<
", "<< stepl/
mat.X0 <<
", "<<
mat.X0 <<
", "<<
mat.L0 <<
", "<<
mat.A <<
", "<<
mat.Z <<
", "<<
mat.rho );
336 m_accumulatedMaterialXX0 += stepl/
mat.X0;
337 m_accumulatedRhoS += stepl*
mat.rho;
338 ATH_MSG_VERBOSE(
"- accumulated material information X/x0 = " << m_accumulatedMaterialXX0);
339 ATH_MSG_VERBOSE(
"- accumulated effective densitity rho*s = " << m_accumulatedRhoS);
342 if (m_mapMaterial && !m_materialMapper.empty()){
344 m_materialMapper->recordMaterialHit(am, positionOnLayer);
349 ATH_MSG_WARNING(
"- associate hit - the layer with index " <<
layer->layerIndex().value() <<
" was not found - should not happen!");
360 if (!propSet)
return;
366 if (confinedLayers) {
374 if (
layer && (*layer).layerIndex().value() ) {
375 ATH_MSG_INFO(
" > LayerIndex: "<< (*layer).layerIndex() );
379 auto curIt = propSet->find((*layer).layerIndex());
380 if (curIt != propSet->end()) {
381 ATH_MSG_INFO(
"LayerMaterial assigned for Layer with index: "<< (*layer).layerIndex() );
383 layer->assignMaterialProperties(*((*curIt).second), 1.);
392 if (confinedVolumes) {
395 ATH_MSG_INFO(
"--> found : "<< volumes.size() <<
"confined TrackingVolumes");
397 for (
const auto & volume : volumes) {
400 assignLayerMaterialProperties(*volume, propSet);
409 ATH_MSG_INFO(
"========================================================================================= ");
412 #ifdef TRKDETDESCR_MEMUSAGE
413 m_memoryLogger.refresh(getpid());
414 ATH_MSG_INFO(
"[ memory usage ] Start building of material maps: " );
419 std::map< std::string, Trk::LayerMaterialMap* > layerMaterialMaps;
420 for (
auto& lmcIter : m_layerMaterialCreators ){
421 ATH_MSG_INFO(
"-> Creating material map '"<< lmcIter->layerMaterialName() <<
"' from creator "<< lmcIter.typeAndName() );
425 ATH_MSG_INFO( m_layerRecords.size() <<
" LayerMaterialRecords to be finalized for this material mapping run.");
428 for (
auto& lIter : m_layerRecords ) {
435 std::string
vName = eVolume ? (eVolume->
volumeName()) :
" BoundaryCollection ";
438 (lIter.second).finalizeRun(m_mapComposition);
440 if (!m_layerMaterialRecordAnalyser.empty() && m_layerMaterialRecordAnalyser->analyseLayerMaterial(*
layer, lIter.second).isFailure() )
443 bool analyse = (m_layerMaterialCreators.size() == m_layerMaterialAnalysers.size());
446 for (
auto& lmcIter : m_layerMaterialCreators ){
448 #ifdef TRKDETDESCR_MEMUSAGE
449 m_memoryLogger.refresh(getpid());
450 ATH_MSG_INFO(
"[ memory usage ] Before building the map for Layer "<< layerKey.
value() );
454 #ifdef TRKDETDESCR_MEMUSAGE
455 m_memoryLogger.refresh(getpid());
456 ATH_MSG_INFO(
"[ memory usage ] After building the map for Layer "<< layerKey.
value() );
462 (*layerMaterialMaps[lmcIter->layerMaterialName()])[layerKey.
value()] = lMaterial;
464 if (analyse && lMaterial && (m_layerMaterialAnalysers[ilmc]->analyseLayerMaterial(*
layer, *lMaterial)).isFailure() )
465 ATH_MSG_WARNING(
"Could not analyse created LayerMaterialProperties for layer "<< layerKey.
value() );
470 ATH_MSG_INFO(
"Finalize map synchronization and write the maps to the DetectorStore.");
472 for (
auto& ilmIter : layerMaterialMaps ){
474 if (m_mapComposition){
476 ilmIter.second->updateElementTable(tElementTable);
477 if (ilmIter.second->elementTable()){
478 ATH_MSG_INFO(
"ElementTable for LayerMaterialMap '" << ilmIter.first <<
"' found and syncrhonized." );
483 if ( (
detStore()->record(ilmIter.second, ilmIter.first,
false)).isFailure()){
484 ATH_MSG_ERROR(
"Writing of LayerMaterialMap with name '" << ilmIter.first <<
"' was not successful." );
485 delete ilmIter.second;
486 }
else ATH_MSG_INFO(
"LayerMaterialMap: " << ilmIter.first <<
" written to the DetectorStore!" );
488 delete m_elementTable;
490 #ifdef TRKDETDESCR_MEMUSAGE
491 m_memoryLogger.refresh(getpid());
492 ATH_MSG_INFO(
"[ memory usage ] At the end of the material map creation.");
496 ATH_MSG_INFO(
"========================================================================================= " );
498 double unmapped = (m_unmapped+m_mapped) ?
double(m_unmapped)/
double(m_unmapped+m_mapped) : 0.;
500 ATH_MSG_INFO(
" -> Skipped (outisde) : " << m_skippedOutside );
501 ATH_MSG_INFO(
"========================================================================================= " );
503 return StatusCode::SUCCESS;
513 m_mappingVolume = trackingGeometry().trackingVolume(m_mappingVolumeName);
516 registerVolume(*trackingVolume, 0);
518 ATH_MSG_INFO(
"Add "<< m_layerRecords.size() <<
" confined volume layers to mapping setup.");
519 ATH_MSG_INFO(
"Add "<< trackingGeometry().numBoundaryLayers() <<
" boundary layers to mapping setup.");
522 for (
const auto bLayerIter : trackingGeometry().boundaryLayers())
523 insertLayerMaterialRecord(*(bLayerIter.first));
525 ATH_MSG_INFO(
"Map for "<< m_layerRecords.size() <<
" layers booked & prepared for mapping procedure");
527 return StatusCode::SUCCESS;
534 int sublevel = lvl+1;
538 std::cout <<
"TrackingVolume name: "<< tvol.
volumeName() << std::endl;
541 std::vector<const Trk::Layer*> volumeLayers;
545 if (confinedLayers) {
550 std::cout <<
"- found : "<<
layers.size() <<
"confined Layers"<< std::endl;
552 auto clIter =
layers.begin();
553 auto clIterE =
layers.end();
554 for ( ; clIter != clIterE; ++clIter ) {
556 const Amg::Vector3D& sReferencePoint = (*clIter)->surfaceRepresentation().globalReferencePoint();
557 bool insideMappingVolume = m_mappingVolume ? m_mappingVolume->inside(sReferencePoint) :
true;
558 if ((*clIter)->layerMaterialProperties() && insideMappingVolume)
559 volumeLayers.push_back((*clIter));
564 for (
auto& lIter : volumeLayers )
565 insertLayerMaterialRecord(*lIter);
569 if (confinedVolumes) {
574 std::cout <<
"- found : "<< volumes.size() <<
"confined TrackingVolumes"<< std::endl;
576 auto volumesIter = volumes.begin();
577 for (; volumesIter != volumes.end(); ++volumesIter)
579 registerVolume(**volumesIter, sublevel);
596 if (layerMaterialBinUtility){
599 layerMaterialBinUtility,
602 m_layerRecords[&lay] = lmr;
607 std::stringstream
msg;
608 msg <<
"Failed to get conditions data " << m_trackingGeometryReadKey.key() <<
".";
609 throw std::runtime_error(
msg.str());