diff --git a/Tracking/TrkDetDescr/TrkDetDescrAlgs/TrkDetDescrAlgs/LayerMaterialRecord.h b/Tracking/TrkDetDescr/TrkDetDescrAlgs/TrkDetDescrAlgs/LayerMaterialRecord.h
deleted file mode 100755
index 8de03f8a0acf95d498deef2944f4c887f76c3737..0000000000000000000000000000000000000000
--- a/Tracking/TrkDetDescr/TrkDetDescrAlgs/TrkDetDescrAlgs/LayerMaterialRecord.h
+++ /dev/null
@@ -1,133 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-//////////////////////////////////////////////////////////////////
-// LayerMaterialRecord.h, (c) ATLAS Detector software
-///////////////////////////////////////////////////////////////////
-
-#ifndef TRKDETDESCRALGS_LAYERMATERIALRECORD_H
-#define TRKDETDESCRALGS_LAYERMATERIALRECORD_H
-
-#include "TrkDetDescrUtils/BinUtility.h"
-#include "TrkDetDescrUtils/MaterialAssociationType.h"
-#include "TrkGeometry/MaterialProperties.h"
-#include "GeoPrimitives/GeoPrimitives.h"
-
-
-namespace Trk {
-
- /** @class LayerMaterialRecord
-
- Helper Class to record the material during the GeantinoNtupleMappingProcess
-
- @author Andreas.Salzburger@cern.ch
- */
-
- class Layer;
- class AssociatedMaterial;
-
- class LayerMaterialRecord {
-
- public :
- /** Default Constructor */
- LayerMaterialRecord();
-
- /** Constructor */
- LayerMaterialRecord(double tlayer,
- const BinUtility* binutils,
- MaterialAssociationType assoc=EffectiveNumAtoms,
- int maxSiliconHits = 0);
-
- /** Constructor */
- LayerMaterialRecord(const LayerMaterialRecord& lmr);
-
- /** Destructor */
- ~LayerMaterialRecord();
-
- /** Assignment operator */
- LayerMaterialRecord& operator=(const LayerMaterialRecord& lmr);
-
- /** adding the information about the Geantino hit */
- void associateGeantinoHit(const Amg::Vector3D& pos,
- double s,
- double x0,
- double l0,
- double a,
- double z,
- double rho,
- int screenoutput=0);
-
- /** finalize the Event */
- AssociatedMaterial* finalizeEvent(const Trk::Layer& lay, int screenoutput=0);
-
- /** finalize the Run */
- void finalizeRun(int screenoutput=0);
-
- /** return method for the LayerMaterial */
- const MaterialPropertiesMatrix& associatedLayerMaterial() const;
-
- /** return the BinUtility */
- const Trk::BinUtility* binUtility() const;
-
- private:
-
- double m_layerThickness; //!< record the layerThickness
- BinUtility* m_binUtility; //!< record the BinnedArray
- int m_bins1; //!< number of bins in coordinate 1
- int m_bins2; //!< number of bins in coordinate 2
-
- //!< event normalizers
- int m_steps;
- double m_track_effNumAtoms;
-
- //!< event related information
- double m_s;
- double m_s_x0;
- double m_s_l0;
- double m_a;
- double m_z;
- double m_rho;
-
- double m_rPosition;
- double m_zPosition;
- double m_phiPosition;
- double m_thetaPosition;
-
- MaterialAssociationType m_assoc; //!< type of hit association
-
- //!< run normalizers
- std::vector< std::vector<int> > m_events;
-
- //!< event related information
- std::vector< std::vector<double> > m_run_t_x0;
- std::vector< std::vector<double> > m_run_t_l0;
- std::vector< std::vector<double> > m_run_a;
- std::vector< std::vector<double> > m_run_z;
- std::vector< std::vector<double> > m_run_rho;
-
- int m_maxSiliconHits;
- int m_siliconHits;
-
- //!< the final material properties
- mutable MaterialPropertiesMatrix m_associatedLayerMaterial;
-
- static double s_siliconRadiationLength;
-
- //!< clear the material -> calls delete
- void clearMaterial();
- //!< copy from another vector
- void copyMaterial(const MaterialPropertiesMatrix& mat);
-
- };
-
- inline const MaterialPropertiesMatrix& LayerMaterialRecord::associatedLayerMaterial() const
- { return m_associatedLayerMaterial; }
-
- inline const Trk::BinUtility* LayerMaterialRecord::binUtility() const
- { return m_binUtility; }
-
-}
-
-#endif
-
diff --git a/Tracking/TrkDetDescr/TrkDetDescrAlgs/TrkDetDescrAlgs/MaterialMapping.h b/Tracking/TrkDetDescr/TrkDetDescrAlgs/TrkDetDescrAlgs/MaterialMapping.h
index d628752fe3beddc3fc6698bf8bb4cd17f87543a5..693bc01cdbca98190958858a31951db7c802399c 100644
--- a/Tracking/TrkDetDescr/TrkDetDescrAlgs/TrkDetDescrAlgs/MaterialMapping.h
+++ b/Tracking/TrkDetDescr/TrkDetDescrAlgs/TrkDetDescrAlgs/MaterialMapping.h
@@ -20,6 +20,20 @@
#include <map>
#include "GeoPrimitives/GeoPrimitives.h"
+#ifdef TRKDETDESCR_MEMUSAGE
+#include "TrkDetDescrUtils/MemoryLogger.h"
+#endif
+
+#ifndef UCHARCONV
+#define UCHARCONV
+#define ucharbin 0.00392157 // 1./double(1.*UCHAR_MAX)
+// int to unsigned char and vv
+#define uchar2uint(uchar) static_cast<unsigned int>(uchar)
+#define uint2uchar(unint) static_cast<unsigned char>(unint)
+// double to unsigned char and vv
+#define uchar2dfrac(uchar) double(uchar*ucharbin)
+#define dfrac2uchar(dfrac) lrint(dfrac*UCHAR_MAX)
+#endif
class TTree;
@@ -31,18 +45,15 @@ namespace Trk {
class SurfaceMaterialRecord;
class LayerMaterialRecord;
class LayerMaterialMap;
+ class Material;
class MaterialProperties;
class BinnedLayerMaterial;
class CompressedLayerMaterial;
- class IMaterialMapper;
+ class ElementTable;
+ class IMaterialMapper;
+ class ILayerMaterialAnalyser;
+ class ILayerMaterialCreator;
class ITrackingGeometrySvc;
-
- struct IndexedMaterial {
- unsigned short int firstBin;
- unsigned short int secondBin;
- const Trk::MaterialProperties* materialProperties;
- };
-
/** @class MaterialMapping
@@ -59,6 +70,7 @@ namespace Trk {
/** Standard Athena-Algorithm Constructor */
MaterialMapping(const std::string& name, ISvcLocator* pSvcLocator);
+
/** Default Destructor */
~MaterialMapping();
@@ -77,11 +89,7 @@ namespace Trk {
bool associateHit(const Trk::TrackingVolume& tvol,
const Amg::Vector3D& pos,
double stepl,
- double x0,
- double l0,
- double a,
- double z,
- double rho);
+ const Trk::Material& mat);
/** Output information with Level */
void registerVolume(const Trk::TrackingVolume& tvol, int lvl);
@@ -89,8 +97,8 @@ namespace Trk {
//!< private void Method to map the layer material
void assignLayerMaterialProperties(const Trk::TrackingVolume& tvol, Trk::LayerMaterialMap* propSet);
- //!< private method to compress material
- Trk::CompressedLayerMaterial* compressMaterial(const Trk::BinnedLayerMaterial& binnedLayerMaterial);
+ //!< create the LayerMaterialRecord */
+ void insertLayerMaterialRecord(const Trk::Layer& lay);
/** Retrieve the TrackingGeometry */
StatusCode retrieveTrackingGeometry();
@@ -108,34 +116,37 @@ namespace Trk {
/** general steering */
double m_etaCutOff;
+ bool m_useLayerThickness; //!< use the actual layer thickness
int m_associationType;
+
+ ToolHandle<ILayerMaterialAnalyser> m_layerMaterialRecordAnalyser; //!< the layer material analyser for the layer material record
+ ToolHandleArray<ILayerMaterialAnalyser> m_layerMaterialAnalysers; //!< the layer material analysers per creator (if wanted)
+ ToolHandleArray<ILayerMaterialCreator> m_layerMaterialCreators; //!< the layer material creators
/** Mapper and Inspector */
bool m_mapMaterial;
- ToolHandle<IMaterialMapper> m_materialMapper; //!< Pointer to an IMaterialMapper algTool
- int m_materialMappingEvents; //!< count the number of validation records to avoid 2G files
- int m_maxMaterialMappingEvents; //!< limit the number of validation records to avoid 2G files
-
+ ToolHandle<IMaterialMapper> m_materialMapper; //!< Pointer to an IMaterialMapper algTool
+ bool m_mapComposition; //!< map the composition of the material
+ double m_minCompositionFraction; //!< minimal fraction to be accounted for the composition recording
+
+
+ Trk::ElementTable* m_elementTable; //!< the accumulated element table
+ std::string m_inputEventElementTable; //!< input event table
+
// the material maps ordered with layer keys
std::map<const Layer*, LayerMaterialRecord > m_layerRecords;
-
- // compression
- bool m_compressMaterial;
- double m_compressedMaterialThickness;
- size_t m_compressedMaterialX0Bins;
- size_t m_compressedMaterialZARhoBins;
- // statistics
+ // statistics for steps
size_t m_mapped;
size_t m_unmapped;
size_t m_skippedOutside;
- // break for validation
- size_t m_mappedEvents;
- size_t m_maxMappedEvents;
-
int m_layerMaterialScreenOutput;
+#ifdef TRKDETDESCR_MEMUSAGE
+ MemoryLogger m_memoryLogger; //!< in case the memory is logged
+#endif
+
};
}
diff --git a/Tracking/TrkDetDescr/TrkDetDescrAlgs/TrkDetDescrAlgs/MaterialValidation.h b/Tracking/TrkDetDescr/TrkDetDescrAlgs/TrkDetDescrAlgs/MaterialValidation.h
index c51e1a2a7a2419edd7788bebbbf5b1786440c025..086f121f96ff81c13cd80c1326aef051fdedcbe0 100755
--- a/Tracking/TrkDetDescr/TrkDetDescrAlgs/TrkDetDescrAlgs/MaterialValidation.h
+++ b/Tracking/TrkDetDescr/TrkDetDescrAlgs/TrkDetDescrAlgs/MaterialValidation.h
@@ -26,6 +26,7 @@ namespace Trk {
class IMaterialMapper;
class TrackingGeometry;
class TrackingVolume;
+ class Surface;
typedef std::pair<Amg::Vector3D, const Trk::TrackingVolume*> PositionAtBoundary;
@@ -34,6 +35,23 @@ namespace Trk {
@author Andreas.Salzburger@cern.ch
*/
+
+ // exit the volume
+ struct VolumeExit {
+ // the triple struct
+ const Trk::TrackingVolume* nVolume;
+ const Trk::Surface* bSurface;
+ Amg::Vector3D vExit;
+
+ VolumeExit( const Trk::TrackingVolume* itv = 0,
+ const Trk::Surface* ibs = 0,
+ Amg::Vector3D iexit = Amg::Vector3D(0.,0.,0.)) :
+ nVolume(itv),
+ bSurface(ibs),
+ vExit(iexit){}
+
+ };
+
class MaterialValidation : public AthAlgorithm {
public:
@@ -70,9 +88,13 @@ namespace Trk {
int m_maxMaterialMappingEvents; //!< limit the number of validation records to avoid 2G files
- Rndm::Numbers* m_flatDist; //!< Random generator for flat distribution
- double m_etaMin;
- double m_etaMax;
+ Rndm::Numbers* m_flatDist; //!< Random generator for flat distribution
+ double m_etaMin; //!< eta boundary
+ double m_etaMax; //!< eta boundary
+ bool m_runNativeNavigation; //!< validate the native TG navigation
+
+ double m_accTinX0; //!< accumulated t in X0
+
};
diff --git a/Tracking/TrkDetDescr/TrkDetDescrAlgs/cmt/requirements b/Tracking/TrkDetDescr/TrkDetDescrAlgs/cmt/requirements
index e463d311fe209ad2048ddee7796cbff720b43b19..a94ac811b86c77bbcb3050f7a99e23d6747fe5e6 100755
--- a/Tracking/TrkDetDescr/TrkDetDescrAlgs/cmt/requirements
+++ b/Tracking/TrkDetDescr/TrkDetDescrAlgs/cmt/requirements
@@ -9,25 +9,24 @@ use AtlasPolicy AtlasPolicy-*
use AthenaBaseComps AthenaBaseComps-* Control
use AtlasROOT AtlasROOT-* External
use GaudiInterface GaudiInterface-* External
-use TrkDetDescrUtils TrkDetDescrUtils-* Tracking/TrkDetDescr
-use TrkGeometry TrkGeometry-* Tracking/TrkDetDescr
use GeoPrimitives GeoPrimitives-* DetectorDescription
+use TrkDetDescrUtils TrkDetDescrUtils-* Tracking/TrkDetDescr
#################################################################
# private use statements
private
use TrkDetDescrInterfaces TrkDetDescrInterfaces-* Tracking/TrkDetDescr
-use TrkSurfaces TrkSurfaces-* Tracking/TrkDetDescr
use TrkVolumes TrkVolumes-* Tracking/TrkDetDescr
-
+use TrkGeometry TrkGeometry-* Tracking/TrkDetDescr
+use TrkNeutralParameters TrkNeutralParameters-* Tracking/TrkEvent
public
library TrkDetDescrAlgs *.cxx components/*.cxx
apply_pattern component_library
# uncomment this line if you do want to build the memory monitoring lines
-# apply_tag use_trkdetdescr_memmon
+apply_tag use_trkdetdescr_memmon
macro_append TrkDetDescrMemMonMacro "-DTRKDETDESCR_MEMUSAGE"
diff --git a/Tracking/TrkDetDescr/TrkDetDescrAlgs/src/LayerMaterialRecord.cxx b/Tracking/TrkDetDescr/TrkDetDescrAlgs/src/LayerMaterialRecord.cxx
deleted file mode 100755
index 9d260bda70c2389a057bf946ad8feec53570baff..0000000000000000000000000000000000000000
--- a/Tracking/TrkDetDescr/TrkDetDescrAlgs/src/LayerMaterialRecord.cxx
+++ /dev/null
@@ -1,425 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-//////////////////////////////////////////////////////////////////
-// LayerMaterialRecord.cxx, (c) ATLAS Detector software
-///////////////////////////////////////////////////////////////////
-
-#include "TrkDetDescrAlgs/LayerMaterialRecord.h"
-#include "TrkGeometry/MaterialProperties.h"
-#include "TrkGeometry/AssociatedMaterial.h"
-#include "TrkGeometry/Layer.h"
-#include "TrkGeometry/CylinderLayer.h"
-#include "TrkGeometry/DiscLayer.h"
-#include "TrkSurfaces/Surface.h"
-#include "TrkSurfaces/SurfaceBounds.h"
-#include "TrkSurfaces/DiscBounds.h"
-#include "TrkSurfaces/CylinderBounds.h"
-
-double Trk::LayerMaterialRecord::s_siliconRadiationLength = 93.6;
-
-Trk::LayerMaterialRecord::LayerMaterialRecord()
-: m_layerThickness(0.),
- m_binUtility(0),
- m_bins1(0),
- m_bins2(0),
- m_steps(0),
- m_track_effNumAtoms(0.),
- m_s(0.),
- m_s_x0(0.),
- m_s_l0(0.),
- m_a(0.),
- m_z(0.),
- m_rho(0.),
- m_rPosition(0.),
- m_zPosition(0.),
- m_phiPosition(0.),
- m_thetaPosition(0.),
- m_assoc(Trk::EffectiveNumAtoms),
- m_maxSiliconHits(0),
- m_siliconHits(0)
-{}
-
-Trk::LayerMaterialRecord::LayerMaterialRecord(double thickness,
- const BinUtility* binutils,
- Trk::MaterialAssociationType assoc,
- int maxSiliconHits)
-: m_layerThickness(thickness),
- m_binUtility(binutils ? binutils->clone() : 0),
- m_bins1(binutils ? (binutils->max(0) + 1) : 1),
- m_bins2(binutils && binutils->dimensions() > 1 ? (binutils->max(1) + 1) : 1),
- m_steps(0),
- m_track_effNumAtoms(0.),
- m_s(0.),
- m_s_x0(0.),
- m_s_l0(0.),
- m_a(0.),
- m_z(0.),
- m_rho(0.),
- m_rPosition(0.),
- m_zPosition(0.),
- m_phiPosition(0.),
- m_thetaPosition(0.),
- m_assoc(assoc),
- m_maxSiliconHits(maxSiliconHits),
- m_siliconHits(0)
-
-{
- for (int ibin=0; ibin<m_bins2; ++ibin) {
- // event
- m_events.push_back(std::vector<int>(m_bins1, 0));
- // run
- m_run_t_x0.push_back(std::vector<double>(m_bins1, 0.));
- m_run_t_l0.push_back(std::vector<double>(m_bins1, 0.));
- m_run_a.push_back(std::vector<double>(m_bins1, 0.));
- m_run_z.push_back(std::vector<double>(m_bins1, 0.));
- m_run_rho.push_back(std::vector<double>(m_bins1, 0.));
- }
-}
-
-
-Trk::LayerMaterialRecord::LayerMaterialRecord(const Trk::LayerMaterialRecord& lmr)
-: m_layerThickness(lmr.m_layerThickness),
- m_binUtility(lmr.m_binUtility ? lmr.m_binUtility->clone() : 0),
- m_bins1(lmr.m_bins1),
- m_bins2(lmr.m_bins2),
- m_steps(lmr.m_steps),
- m_track_effNumAtoms(lmr.m_track_effNumAtoms),
- m_s(lmr.m_s),
- m_s_x0(lmr.m_s_x0),
- m_s_l0(lmr.m_s_l0),
- m_a(lmr.m_a),
- m_z(lmr.m_z),
- m_rho(lmr.m_rho),
- m_rPosition(lmr.m_rPosition),
- m_zPosition(lmr.m_zPosition),
- m_phiPosition(lmr.m_phiPosition),
- m_thetaPosition(lmr.m_thetaPosition),
- m_assoc(lmr.m_assoc),
- m_events(lmr.m_events),
- m_run_t_x0(lmr.m_run_t_x0),
- m_run_t_l0(lmr.m_run_t_l0),
- m_run_a(lmr.m_run_a),
- m_run_z(lmr.m_run_z),
- m_run_rho(lmr.m_run_rho),
- m_maxSiliconHits(lmr.m_maxSiliconHits),
- m_siliconHits(lmr.m_siliconHits)
-{
- copyMaterial(lmr.m_associatedLayerMaterial);
-}
-
-
-Trk::LayerMaterialRecord& Trk::LayerMaterialRecord::operator=(const Trk::LayerMaterialRecord& lmr)
-{
- if (this != &lmr) {
- m_layerThickness = lmr.m_layerThickness;
- m_binUtility = lmr.m_binUtility ? lmr.m_binUtility->clone() : 0;
- m_bins1 = lmr.m_bins1;
- m_bins2 = lmr.m_bins2;
- m_assoc = lmr.m_assoc;
- m_steps = lmr.m_steps;
- m_s = lmr.m_s;
- m_s_x0 = lmr.m_s_x0;
- m_s_l0 = lmr.m_s_l0;
- m_a = lmr.m_a;
- m_z = lmr.m_z;
- m_rho = lmr.m_rho;
- m_rPosition = lmr.m_rPosition;
- m_zPosition = lmr.m_zPosition;
- m_phiPosition = lmr.m_phiPosition;
- m_thetaPosition = lmr.m_thetaPosition;
- m_events = lmr.m_events;
- m_track_effNumAtoms = lmr.m_track_effNumAtoms;
- m_run_t_x0 = lmr.m_run_t_x0;
- m_run_t_l0 = lmr.m_run_t_l0;
- m_run_a = lmr.m_run_a;
- m_run_z = lmr.m_run_z;
- m_run_rho = lmr.m_run_rho;
- m_maxSiliconHits = lmr.m_maxSiliconHits;
- m_siliconHits = lmr.m_siliconHits;
- // deal with the material
- clearMaterial();
- copyMaterial( lmr.m_associatedLayerMaterial );
- }
-
- return (*this);
-}
-
-
-Trk::LayerMaterialRecord::~LayerMaterialRecord()
-{
- // don't delete the material -> its given to the outside world
- delete m_binUtility;
-}
-
-void Trk::LayerMaterialRecord::associateGeantinoHit(const Amg::Vector3D& pos,
- double s,
- double x0,
- double l0,
- double a,
- double z,
- double rho,
- int)
-{
-
- // spatial information
- m_rPosition += pos.perp();
- m_zPosition += pos.z();
- m_phiPosition += pos.phi();
- m_thetaPosition += pos.theta();
-
- m_track_effNumAtoms += s*rho/a;
- m_steps++;
- m_s += s;
-
- // must be true for all association methods
- m_s_x0 += s/x0;
- m_s_l0 += s/l0;
- m_rho += rho * s;
-
- if (m_maxSiliconHits && fabs(x0-s_siliconRadiationLength)<0.5) ++m_siliconHits;
-
- if (m_assoc==Trk::EffectiveNumSteps) {
- m_a += a;
- m_z += z;
- } else if (m_assoc==Trk::EffectiveNumAtoms) {
- m_a += s*rho;
- m_z += s*rho * z/a;
- } else if (m_assoc==Trk::RadiationLength) {
- m_a += a*s/x0;
- m_z += z*s/x0;
- } else {
- m_a += a*s;
- m_z += z*s;
- }
-
-}
-
-Trk::AssociatedMaterial* Trk::LayerMaterialRecord::finalizeEvent(const Trk::Layer& lay, int screenoutput)
-{
- // the corrFactor
- double corrFactorInv = 0.;
- // the return value
- Trk::AssociatedMaterial* returnStep = 0;
-
- if (screenoutput)
- std::cout << " o Finalize Event for layer " << lay.layerIndex().value() << std::endl;
-
-
- // the silicon double counting can veto the entire event
- bool vetoEvent = (m_maxSiliconHits && m_siliconHits > m_maxSiliconHits) ? true : false;
-
- if (!vetoEvent && m_s_x0 >10e-10 && m_bins1 && m_bins2 && m_steps) {
-
- // the densed phi/theta values
- double hitPhi = m_phiPosition/m_steps;
- double hitTheta = m_thetaPosition/m_steps;
-
- // the densed hit information
- double densedHitX = 0.;
- double densedHitY = 0.;
- double densedHitZ = 0.;
- double densedHitR = 0.;
-
- // cast it to cylinder/disc
- const Trk::CylinderLayer* cylLayer = dynamic_cast<const Trk::CylinderLayer*>(&lay);
-
- // get the position
- if (cylLayer) {
- // Cylinder Bounds needed
- const Trk::CylinderBounds& cylBounds = cylLayer->surfaceRepresentation().bounds();
- densedHitR = cylBounds.r();
- double halflengthZ = cylBounds.halflengthZ();
- double centerZ = lay.surfaceRepresentation().center().z();
- densedHitZ = densedHitR/tan(hitTheta);
- // the corrfactor is
- corrFactorInv = fabs(sin(hitTheta));
- // veto the event if the densed radius is outside the cylinder
- vetoEvent = densedHitZ < centerZ-halflengthZ || densedHitZ > centerZ+halflengthZ;
-
- } else {
- // DiscBounds needed
- const Trk::DiscBounds* discBounds = dynamic_cast<const Trk::DiscBounds*>(&(lay.surfaceRepresentation().bounds()));
- if (discBounds) {
- densedHitZ = lay.surfaceRepresentation().center().z();
- densedHitR = densedHitZ*tan(hitTheta);
- // the corrfactor is
- corrFactorInv = fabs(cos(hitTheta));
- // veto the event if the densed radius is outside the disc
- vetoEvent = densedHitR < discBounds->rMin() || densedHitR > discBounds->rMax();
- }
- }
-
- if (!vetoEvent) {
-
- densedHitX = densedHitR*cos(hitPhi);
- densedHitY = densedHitR*sin(hitPhi);
-
- Amg::Vector3D hitPosition(densedHitX,densedHitY,densedHitZ);
-
- int recordBin1 = m_binUtility ? m_binUtility->bin(hitPosition, 0) : 0;
- int recordBin2 = ( m_binUtility && m_binUtility->dimensions() > 1) ? m_binUtility->bin(hitPosition, 1) : 0;
-
- // safety checks
- recordBin1 = ( recordBin1 < 0 ) ? 0 : recordBin1;
- recordBin1 = ( recordBin1 >= m_bins1) ? m_bins1-1 : recordBin1;
- recordBin2 = ( recordBin2 < 0 ) ? 0 : recordBin2;
- recordBin2 = ( recordBin2 >= m_bins2) ? m_bins2-1 : recordBin2;
-
- // finalization
- m_rho *= corrFactorInv;
- m_rho /= (m_layerThickness); // condense to layer thickness
- // the different methods
- if (m_assoc==Trk::EffectiveNumSteps && m_steps) {
- m_a /= m_steps;
- m_z /= m_steps;
- } else if (m_assoc==Trk::EffectiveNumAtoms && m_track_effNumAtoms>10e-23) {
- m_a /= m_track_effNumAtoms;
- m_z /= m_track_effNumAtoms;
- } else if (m_assoc==Trk::RadiationLength && m_s_x0 > 10e-10) {
- m_a /= m_s_x0;
- m_z /= m_s_x0;
- } else if (m_s > 10e-10) {
- m_a /= m_s;
- m_z /= m_s;
- }
-
- // assign the variables for the outside world
- double x0 = m_s/m_s_x0;
- double l0 = m_s/m_s_l0;
-
- if (screenoutput > 1)
- std::cout << " Collected information for bin ["
- << recordBin1 << ", " << recordBin2 << "] in this event ( steps, s/x0, x0 ) = "
- << m_steps << ", " << m_s_x0 << ", " << x0 << std::endl;
-
- // assign the run variables
- m_events[recordBin2][recordBin1]++;
-
- m_run_t_x0[recordBin2][recordBin1] += m_s_x0*corrFactorInv;
- m_run_t_l0[recordBin2][recordBin1] += m_s_l0*corrFactorInv;
- m_run_a[recordBin2][recordBin1] += m_a;
- m_run_z[recordBin2][recordBin1] += m_z;
- m_run_rho[recordBin2][recordBin1] += m_rho;
-
- if (screenoutput >1)
- std::cout << " Cummulated material from " << m_events[recordBin2][recordBin1]
- << " events ( s/x0 ) = " << m_run_t_x0[recordBin2][recordBin1] << std::endl;
-
- // create the return step
- returnStep = (corrFactorInv >1.) ? 0 : new Trk::AssociatedMaterial( hitPosition,
- m_s_x0*x0,
- x0,
- l0,
- m_a,
- m_z,
- m_rho,
- 1./corrFactorInv,
- lay.enclosingTrackingVolume(),
- &lay);
-
- }
- }
-
- // reset silicon hits
- m_siliconHits = 0;
-
- // reset event variables
- m_steps = 0;
- m_s = 0.;
- m_s_x0 = 0.;
- m_s_l0 = 0.;
- m_a = 0.;
- m_z = 0.;
- m_rho = 0.;
- m_track_effNumAtoms = 0.;
-
- m_rPosition = 0.;
- m_zPosition = 0.;
- m_phiPosition = 0.;
- m_thetaPosition = 0.;
-
- return returnStep;
-
-}
-
-void Trk::LayerMaterialRecord::finalizeRun(int screenoutput)
-{
-
- if (screenoutput)
- std::cout << " o Finalize Layer with [bins1,bins2] = [" << m_bins1 << ", " << m_bins2 << "] material bins " << std::endl;
-
- m_associatedLayerMaterial.reserve(m_bins2);
-
- for (int ibin2 = 0; ibin2 < m_bins2; ++ibin2) {
- // create the vector first
- Trk::MaterialPropertiesVector matVector;
- matVector.reserve(m_bins1);
- // loop over local 1 bins
- for (int ibin1 = 0; ibin1 < m_bins1; ++ibin1) {
- Trk::MaterialProperties* binMaterial = 0;
- if (m_events[ibin2][ibin1]) {
-
- if (screenoutput)
- std::cout << " Bin [" << ibin1 << "," << ibin2 << "]" << " has " << m_events[ibin2][ibin1] << " recorded events." << std::endl;
-
- m_run_t_x0[ibin2][ibin1] /= m_events[ibin2][ibin1];
- m_run_t_l0[ibin2][ibin1] /= m_events[ibin2][ibin1];
- m_run_a[ibin2][ibin1] /= m_events[ibin2][ibin1];
- m_run_z[ibin2][ibin1] /= m_events[ibin2][ibin1];
- m_run_rho[ibin2][ibin1] /= m_events[ibin2][ibin1];
-
- std::cout << " Averaged material from " << m_events[ibin2][ibin1]
- << " events ( s/x0 ) = " << m_run_t_x0[ibin2][ibin1] << std::endl;
-
- binMaterial = new Trk::MaterialProperties(1.,
- 1./m_run_t_x0[ibin2][ibin1],
- 1./m_run_t_l0[ibin2][ibin1],
- m_run_a[ibin2][ibin1],
- m_run_z[ibin2][ibin1],
- m_run_rho[ibin2][ibin1]);
- if (screenoutput > 1)
- std::cout << " Created Material : " << *binMaterial << std::endl;
-
- }
- matVector.push_back(binMaterial);
- }
- m_associatedLayerMaterial.push_back(matVector);
- }
-}
-
-void Trk::LayerMaterialRecord::clearMaterial()
-{
-
- Trk::MaterialPropertiesMatrix::iterator matMatrixIter = m_associatedLayerMaterial.begin();
- Trk::MaterialPropertiesMatrix::iterator matMatrixIterEnd = m_associatedLayerMaterial.end();
- for ( ; matMatrixIter != matMatrixIterEnd; ++matMatrixIter) {
- // loop over the subsets
- std::vector<const Trk::MaterialProperties*>::iterator matIter = (*matMatrixIter).begin();
- std::vector<const Trk::MaterialProperties*>::iterator matIterEnd = (*matMatrixIter).end();
- for ( ; matIter != matIterEnd; ++matIter)
- delete (*matIter);
- }
-}
-
-void Trk::LayerMaterialRecord::copyMaterial(const MaterialPropertiesMatrix& materialMatrix)
-{
- // clear the vector
- m_associatedLayerMaterial.clear();
-
- Trk::MaterialPropertiesMatrix::const_iterator matMatrixIter = materialMatrix.begin();
- Trk::MaterialPropertiesMatrix::const_iterator matMatrixIterEnd = materialMatrix.end();
- for ( ; matMatrixIter != matMatrixIterEnd; ++matMatrixIter) {
- Trk::MaterialPropertiesVector matProp;
- // loop over the subsets
- std::vector<const Trk::MaterialProperties*>::const_iterator matIter = (*matMatrixIter).begin();
- std::vector<const Trk::MaterialProperties*>::const_iterator matIterEnd = (*matMatrixIter).end();
- for ( ; matIter != matIterEnd; ++matIter) {
- // test it
- matProp.push_back( ( (*matIter) ? (*matIter)->clone() : 0 ) );
- }
- // and now push back the vector
- m_associatedLayerMaterial.push_back(matProp);
- }
-}
diff --git a/Tracking/TrkDetDescr/TrkDetDescrAlgs/src/MaterialMapping.cxx b/Tracking/TrkDetDescr/TrkDetDescrAlgs/src/MaterialMapping.cxx
index c731dafa59c3e86da00f9d4e39bb664bced025ef..d927cb4edb27d63f161fcb3e70ea8cf312ff0b06 100755
--- a/Tracking/TrkDetDescr/TrkDetDescrAlgs/src/MaterialMapping.cxx
+++ b/Tracking/TrkDetDescr/TrkDetDescrAlgs/src/MaterialMapping.cxx
@@ -8,16 +8,18 @@
// Gaudi Units
#include "GaudiKernel/SystemOfUnits.h"
-//TrkDetDescrAlgs
+//TrkDetDescr Algs, Interfaces, Utils
#include "TrkDetDescrAlgs/MaterialMapping.h"
-#include "TrkDetDescrAlgs/LayerMaterialRecord.h"
-// Trk
#include "TrkDetDescrInterfaces/ITrackingGeometrySvc.h"
#include "TrkDetDescrInterfaces/IMaterialMapper.h"
+#include "TrkDetDescrInterfaces/ILayerMaterialCreator.h"
+#include "TrkDetDescrInterfaces/ILayerMaterialAnalyser.h"
#include "TrkDetDescrUtils/GeometryStatics.h"
#include "TrkDetDescrUtils/LayerIndex.h"
#include "TrkDetDescrUtils/BinUtility.h"
-#include "TrkGeometry/EntryLayerProvider.h"
+// TrkGeometry
+#include "TrkGeometry/LayerMaterialRecord.h"
+//#include "TrkGeometry/ElementTable.h"
#include "TrkGeometry/TrackingGeometry.h"
#include "TrkGeometry/TrackingVolume.h"
#include "TrkGeometry/MaterialStep.h"
@@ -34,7 +36,6 @@
#include <unistd.h>
#endif
-
Trk::MaterialMapping::MaterialMapping(const std::string& name, ISvcLocator* pSvcLocator)
: AthAlgorithm(name,pSvcLocator),
m_trackingGeometrySvc("AtlasTrackingGeometrySvc",name),
@@ -42,42 +43,48 @@ Trk::MaterialMapping::MaterialMapping(const std::string& name, ISvcLocator* pSvc
m_mappingVolumeName("Atlas"),
m_mappingVolume(0),
m_inputMaterialStepCollection("MaterialStepRecords"),
- m_outputLayerMaterialSetName("AtlasTrackingMaterial"),
m_etaCutOff(6.0),
+ m_useLayerThickness(false),
m_associationType(1),
+ m_layerMaterialRecordAnalyser(""),
m_mapMaterial(true),
- m_materialMapper("Trk::MaterialMapper/MappingMaterialMapper"),
- m_materialMappingEvents(0),
- m_maxMaterialMappingEvents(25000),
+ m_materialMapper(""),
+ m_mapComposition(true),
+ m_minCompositionFraction(0.005),
+ m_elementTable(0),
+ m_inputEventElementTable("ElementTable"),
m_mapped(0),
m_unmapped(0),
m_skippedOutside(0),
- m_mappedEvents(0),
- m_maxMappedEvents(-1),
m_layerMaterialScreenOutput(0)
+#ifdef TRKDETDESCR_MEMUSAGE
+ ,m_memoryLogger()
+#endif
{
// the name of the TrackingGeometry to be retrieved
declareProperty("TrackingGeometrySvc" , m_trackingGeometrySvc);
declareProperty("MappingVolumeName" , m_mappingVolumeName);
// general steering
declareProperty("EtaCutOff" , m_etaCutOff);
+ // for the analysis of the material
+ declareProperty("LayerMaterialRecordAnalyser" , m_layerMaterialRecordAnalyser);
+ // for the creation of the material
+ declareProperty("LayerMaterialCreators" , m_layerMaterialCreators);
+ declareProperty("LayerMaterialAnalysers" , m_layerMaterialAnalysers);
// the toolhandle of the MaterialMapper to be used
- declareProperty("MapMaterial" , m_mapMaterial);
- declareProperty("MaterialMapper" , m_materialMapper);
- declareProperty("MaximumMappingEvents" , m_maxMaterialMappingEvents);
- // Compression section
- declareProperty("CompressMaterial" , m_compressMaterial);
- declareProperty("CompressMaterialThickness" , m_compressedMaterialThickness);
- declareProperty("CompressMaterialBinsX0" , m_compressedMaterialX0Bins);
- declareProperty("CompressMaterialBinsZARho" , m_compressedMaterialZARhoBins);
+ declareProperty("MapMaterial" , m_mapMaterial);
+ declareProperty("MaterialMapper" , m_materialMapper);
+ // Composition related parameters
+ declareProperty("MapComposition" , m_mapComposition);
+ declareProperty("MinCompositionFraction" , m_minCompositionFraction);
+ // Steer the layer thickness
+ declareProperty("UseActualLayerThicknesss" , m_useLayerThickness);
// some job setup
- declareProperty("MaterialAssociationType" , m_associationType);
- declareProperty("InputMaterialStepCollection" , m_inputMaterialStepCollection);
- declareProperty("OutputLayerMaterialSetName" , m_outputLayerMaterialSetName);
- // Maximum of mapped events
- declareProperty("MaximumMappedEvents" , m_maxMappedEvents);
+ declareProperty("MaterialAssociationType" , m_associationType);
+ declareProperty("InputMaterialStepCollection" , m_inputMaterialStepCollection);
+ declareProperty("InputElementTable" , m_inputEventElementTable);
// Output screen stuff
- declareProperty("MaterialScreenOutputLevel" , m_layerMaterialScreenOutput);
+ declareProperty("MaterialScreenOutputLevel" , m_layerMaterialScreenOutput);
}
@@ -88,13 +95,22 @@ Trk::MaterialMapping::~MaterialMapping()
StatusCode Trk::MaterialMapping::initialize()
{
- ATH_MSG_INFO( "initialize()" );
+ ATH_MSG_INFO("initialize()");
if ( (m_trackingGeometrySvc.retrieve()).isFailure() )
ATH_MSG_WARNING("Could not retrieve TrackingGeometrySvc");
- if ( (m_materialMapper.retrieve()).isFailure() )
+ if (!m_materialMapper.empty() && (m_materialMapper.retrieve()).isFailure() )
ATH_MSG_WARNING("Could not retrieve MaterialMapper");
+
+ if ( !m_layerMaterialRecordAnalyser.empty() && m_layerMaterialRecordAnalyser.retrieve().isFailure() )
+ ATH_MSG_WARNING("Could not retrieve LayerMaterialAnalyser");
+
+ if ( m_layerMaterialCreators.size() && m_layerMaterialCreators.retrieve().isFailure() )
+ ATH_MSG_WARNING("Could not retrieve any LayerMaterialCreators");
+
+ if ( m_layerMaterialAnalysers.size() && m_layerMaterialAnalysers.retrieve().isFailure() )
+ ATH_MSG_WARNING("Could not retrieve any LayerMaterialAnalysers");
return StatusCode::SUCCESS;
}
@@ -102,122 +118,79 @@ StatusCode Trk::MaterialMapping::initialize()
StatusCode Trk::MaterialMapping::execute()
{
- ATH_MSG_VERBOSE( "MaterialMapping execute() start" );
+ ATH_MSG_VERBOSE("MaterialMapping execute() start");
// ------------------------------- get the trackingGeometry at first place
if (!m_trackingGeometry) {
StatusCode retrieveCode = retrieveTrackingGeometry();
if (retrieveCode.isFailure()){
- ATH_MSG_INFO( "Could not retrieve TrackingGeometry. Exiting." );
+ ATH_MSG_INFO("Could not retrieve TrackingGeometry. Exiting.");
return retrieveCode;
}
- }
+ }
+ if (m_trackingGeometry)
+ ATH_MSG_VERBOSE("TrackingGeometry sucessfully retrieved");
// get the material step vector
const MaterialStepCollection* materialSteps = 0;
-
- // try this
StatusCode retrieveCode = evtStore()->retrieve(materialSteps,m_inputMaterialStepCollection);
if (retrieveCode.isFailure()) {
-
- ATH_MSG_WARNING( "[!] Could not retrieve the step vector!");
+ ATH_MSG_WARNING("[!] Could not retrieve the step vector!");
return retrieveCode;
-
- } else if (m_maxMappedEvents > 0 && m_mappedEvents >= m_maxMappedEvents ) {
-
- ATH_MSG_VERBOSE( "Maximunm mapped events already reached. Skipping.");
- return StatusCode::SUCCESS;
-
} else { // successful retrieval
- // average numbers for the weighting
- int associatedSteps = 0 ;
- double av_eta = 0.;
- double av_phi = 0.;
-
- double eta = 0.;
- double phi = 0.;
-
- ATH_MSG_VERBOSE( "[+] Successfully read " << materialSteps->size() << " geantino steps" );
-
- auto stepIter = materialSteps->begin();
- auto stepIterEnd = materialSteps->end();
+ const Trk::ElementTable* eTableEvent = 0;
+ if (evtStore()->contains<Trk::ElementTable>(m_inputEventElementTable) && evtStore()->retrieve(eTableEvent,m_inputEventElementTable).isFailure()){
+ ATH_MSG_WARNING("No ElementTable coud be retrieved. Switching composition recording off.");
+ } else if (eTableEvent)
+ (*m_elementTable) += (*eTableEvent); // accummulate the table
+ m_mapComposition = eTableEvent ? true : false;
+
+ ATH_MSG_DEBUG("[+] Successfully read "<< materialSteps->size() << " geantino steps");
+ int associatedSteps = 0;
// loop over the hits and associate that to the layers
- for ( ; stepIter != stepIterEnd; ++stepIter) {
-
- double hitX = (*stepIter)->hitX();
- double hitY = (*stepIter)->hitY();
- double hitZ = (*stepIter)->hitZ();
-
- double t = (*stepIter)->steplength();
- double x0 = (*stepIter)->x0();
- double l0 = (*stepIter)->l0();
- double a = (*stepIter)->A();
- double z = (*stepIter)->Z();
- double rho = (*stepIter)->rho();
-
- Amg::Vector3D pos(hitX,hitY,hitZ);
-
- // skip if outside the mapping volume
- if (m_mappingVolume && !m_mappingVolume->inside(pos)){
+ for ( auto& step : (*materialSteps) ) {
+
+ // step length and position
+ double t = step->steplength();
+ Amg::Vector3D pos(step->hitX(), step->hitY(), step->hitZ());
+ // skip if :
+ // -- 0) no mapping volume exists
+ // -- 1) outside the mapping volume
+ // -- 2) outside the eta acceptance
+ if (!m_mappingVolume || !(m_mappingVolume->inside(pos)) || fabs(pos.eta()) > m_etaCutOff ){
++m_skippedOutside;
continue;
}
- // go ahead
- eta = pos.eta();
- phi = pos.phi();
-
// get the lowest TrackingVolume
const Trk::TrackingVolume* associatedVolume = m_trackingGeometry->lowestTrackingVolume(pos);
// add associated volume to volumes map
- if (!associatedVolume)
+ if (!associatedVolume) {
+ ATH_MSG_FATAL("No associated TrackingVolume found. Aborting.");
return StatusCode::FAILURE;
-
- // ------ CUT ------ (e.g. dummy/giant steps, non-physical atomic numbers, high eta ...)
- if (pos.mag() > 1. * Gaudi::Units::mm && t > 10e-10 && t < 10000. && a < 252. && a > 0. && fabs(eta) < m_etaCutOff) {
- if (associateHit(*associatedVolume, pos,t,x0,l0,a,z,rho)) {
- // average over eta, phi, theta
- av_eta += eta;
- av_phi += phi;
- ++associatedSteps;
- }
}
+ // associate the hit
+ if (associateHit(*associatedVolume, pos, t, step->fullMaterial()))
+ ++associatedSteps;
+
} // end of loop over steps
- ATH_MSG_VERBOSE(" Loop over Material steps done, now finalizing event ...");
+ ATH_MSG_VERBOSE("Loop over Material steps done, now finalizing event ...");
// check whether the event was good for at least one hit
if (associatedSteps) {
-
- // at least one step was associated -> thus the event counts
- ++m_mappedEvents;
-
- ATH_MSG_VERBOSE(" There are associated steps, need to call finalizeEvent() & record to the MaterialMapper.");
- // accepted hits
- av_eta /= associatedSteps;
- av_phi /= associatedSteps;
-
+ ATH_MSG_VERBOSE("There are associated steps, need to call finalizeEvent() & record to the MaterialMapper.");
// finalize the event --------------------- Layers ---------------------------------------------
- auto lIter = m_layerRecords.begin();
- for ( ; lIter != m_layerRecords.end(); ++lIter) {
-
- const Trk::Layer* layer = lIter->first;
- if (!layer) continue;
-
- Trk::AssociatedMaterial* assMatHit = (lIter->second).finalizeEvent(*layer, m_layerMaterialScreenOutput);
- if (assMatHit)
- m_materialMapper->recordFullLayerHit(*assMatHit);
+ for (auto& lRecord : m_layerRecords ) {
+ // associated material
+ Trk::AssociatedMaterial* assMatHit = lRecord.second.finalizeEvent((*lRecord.first));
+ // record the full layer hit
+ if (assMatHit && !m_materialMapper.empty()) m_materialMapper->recordLayerHit(*assMatHit, true);
delete assMatHit;
-
// call the material mapper finalize method
- ATH_MSG_VERBOSE( " Calling finalizeEvent on the MaterialMapper ..." );
- // increase the number
- ++m_materialMappingEvents;
- // material mapper
- m_materialMapper->finalizeEvent(av_eta,av_phi);
-
+ ATH_MSG_VERBOSE("Calling finalizeEvent on the MaterialMapper ...");
}
} // the event had at least one associated hit
}
@@ -229,61 +202,31 @@ StatusCode Trk::MaterialMapping::execute()
bool Trk::MaterialMapping::associateHit( const Trk::TrackingVolume& associatedVolume,
const Amg::Vector3D& pos,
double stepl,
- double x0,
- double l0,
- double a,
- double z,
- double rho)
-{
-
-
-
+ const Trk::Material& mat)
+{
// get the intersection with the layer
Amg::Vector3D dir = pos.normalized();
- const Trk::LayerIntersection lIntersect = associatedVolume.closestMaterialLayer(pos, dir, Trk::mappingMode, true);
- if (lIntersect.sIntersection.valid){
+ const Trk::LayerIntersection<Amg::Vector3D> lIntersect = associatedVolume.closestMaterialLayer(pos, dir, Trk::mappingMode, true);
+ if (lIntersect.intersection.valid){
++m_mapped;
// try to find the layer material record
- auto clIter = m_layerRecords.find(lIntersect.layer);
+ auto clIter = m_layerRecords.find(lIntersect.object);
if (clIter != m_layerRecords.end() ){
// record the plain information w/o correction & intersection
- if (m_mapMaterial)
- m_materialMapper->recordMaterialHit(Trk::AssociatedMaterial(pos,
- stepl,
- x0,
- l0,
- a,
- z,
- rho,
- 1.,
- &associatedVolume,
- lIntersect.layer),
- lIntersect.sIntersection.intersection);
+ if (m_mapMaterial && !m_materialMapper.empty())
+ m_materialMapper->recordMaterialHit(Trk::AssociatedMaterial(pos, stepl, mat, 1., &associatedVolume, lIntersect.object), lIntersect.intersection.position);
// LayerMaterialRecord found, add the hit
- (*clIter).second.associateGeantinoHit(lIntersect.sIntersection.intersection, stepl, x0, l0, a, z, rho, m_layerMaterialScreenOutput);
- ATH_MSG_VERBOSE(" - associate Geantino Information ( s/x0 , x0 , l0, a, z, rho ) = "
- << stepl/x0 << ", " << x0 << " , " << l0 << ", " << a << ", " << z << ", " << rho );
- ATH_MSG_VERBOSE(" - to layer with index " << lIntersect.layer->layerIndex().value() << " in '" << associatedVolume.volumeName() << "'.");
+ (*clIter).second.associateGeantinoHit(lIntersect.intersection.position, stepl, mat);
+ ATH_MSG_VERBOSE("- associate Geantino Information at intersection [r/z] = "<< lIntersect.intersection.position.perp() << "/"<< lIntersect.intersection.position.z() );
+ ATH_MSG_VERBOSE("- associate Geantino Information ( s, s/x0 , x0 , l0, a, z, rho ) = "
+ << stepl << ", "<< stepl/mat.X0 << ", "<< mat.X0 << ", "<< mat.L0 << ", "<< mat.A << ", "<< mat.Z << ", "<< mat.rho );
+ ATH_MSG_VERBOSE("- to layer with index "<< lIntersect.object->layerIndex().value() << "in '"<< associatedVolume.volumeName() << "'.");
} else {
ATH_MSG_WARNING("No Layer found in the associated map! Should not happen.");
return false;
}
- } else {
- // unmapped - record for validation
+ } else
++m_unmapped;
- if (m_mapMaterial)
- m_materialMapper->recordMaterialHit(Trk::AssociatedMaterial(pos,
- stepl,
- x0,
- l0,
- a,
- z,
- rho,
- 1.,
- &associatedVolume,
- 0), Amg::Vector3D(0.,0.,0.));
-
- }
// return
return true;
}
@@ -295,60 +238,28 @@ void Trk::MaterialMapping::assignLayerMaterialProperties( const Trk::TrackingVol
if (!propSet) return;
- ATH_MSG_INFO( "Processing TrackingVolume: " << tvol.volumeName() );
-
- // ----------------------------------- loop over the entry layers ----------------------------
- const Trk::EntryLayerProvider* entryLayerProvider = tvol.entryLayerProvider();
- if (entryLayerProvider) {
- // get the objects in a vector-like format
- const std::vector<const Trk::Layer*>& elayers = entryLayerProvider->layers();
- ATH_MSG_INFO( " --> found : " << elayers.size() << " entry Layers" );
- // iterator
- auto elayerIter = elayers.begin();
- // loop over layers
- for ( ; elayerIter != elayers.end(); ++elayerIter) {
- // check if the layer has material properties
- if (*elayerIter) {
- ATH_MSG_INFO( " > LayerIndex: " << (**elayerIter).layerIndex() );
- // set the material!
- if (propSet) {
- // find the pair
- auto curIt = propSet->find((**elayerIter).layerIndex());
- if (curIt != propSet->end()) {
- ATH_MSG_INFO( "LayerMaterial assigned for Layer with index: " << (**elayerIter).layerIndex() );
- // set it to the layer
- (**elayerIter).assignMaterialProperties(*((*curIt).second), 1.);
-
- }
- }
- }
- }
- }
- // -------------------------------------- END OF ENTRY LAYER SECTION --------------------------------------
+ ATH_MSG_INFO("Processing TrackingVolume: "<< tvol.volumeName() );
-
- // ----------------------------------- loop over confined layers -----------------------------
+ // ----------------------------------- loop over confined layers ------------------------------------------
const Trk::BinnedArray< Trk::Layer >* confinedLayers = tvol.confinedLayers();
if (confinedLayers) {
// get the objects in a vector-like format
const std::vector<const Trk::Layer*>& layers = confinedLayers->arrayObjects();
- ATH_MSG_INFO( " --> found : " << layers.size() << " confined Layers" );
+ ATH_MSG_INFO("--> found : "<< layers.size() << "confined Layers");
// the iterator over the vector
- auto layerIter = layers.begin();
// loop over layers
- for ( ; layerIter != layers.end(); ++layerIter) {
+ for (auto& layer : layers) {
// assign the material and output
- if (*layerIter && (**layerIter).layerIndex().value() ) {
- ATH_MSG_INFO( " > LayerIndex: " << (**layerIter).layerIndex() );
+ if (layer && (*layer).layerIndex().value() ) {
+ ATH_MSG_INFO(" > LayerIndex: "<< (*layer).layerIndex() );
// set the material!
if (propSet) {
// find the pair
- auto curIt = propSet->find((**layerIter).layerIndex());
+ auto curIt = propSet->find((*layer).layerIndex());
if (curIt != propSet->end()) {
- ATH_MSG_INFO( "LayerMaterial assigned for Layer with index: " << (**layerIter).layerIndex() );
+ ATH_MSG_INFO("LayerMaterial assigned for Layer with index: "<< (*layer).layerIndex() );
// set it to the layer
- (**layerIter).assignMaterialProperties(*((*curIt).second), 1.);
-
+ (*layer).assignMaterialProperties(*((*curIt).second), 1.);
}
}
}
@@ -360,14 +271,12 @@ void Trk::MaterialMapping::assignLayerMaterialProperties( const Trk::TrackingVol
if (confinedVolumes) {
// get the objects in a vector-like format
const std::vector<const Trk::TrackingVolume*>& volumes = confinedVolumes->arrayObjects();
- ATH_MSG_INFO( " --> found : " << volumes.size() << " confined TrackingVolumes" );
- // the iterator over the vector
- auto volumeIter = volumes.begin();
+ ATH_MSG_INFO("--> found : "<< volumes.size() << "confined TrackingVolumes");
// loop over volumes
- for ( ; volumeIter != volumes.end(); ++volumeIter) {
+ for (auto& volume : volumes) {
// assing the material and output
- if (*volumeIter)
- assignLayerMaterialProperties(**volumeIter, propSet); // call itself recursively
+ if (volume)
+ assignLayerMaterialProperties(*volume, propSet); // call itself recursively
}
}
}
@@ -376,76 +285,101 @@ void Trk::MaterialMapping::assignLayerMaterialProperties( const Trk::TrackingVol
StatusCode Trk::MaterialMapping::finalize()
{
- ATH_MSG_INFO( "========================================================================================= " );
- ATH_MSG_INFO( "finalize() starts ...");
-
- // record the final map to be written to POOL
- Trk::LayerMaterialMap* associatedLayerMaterialProperties = new Trk::LayerMaterialMap();
-
- ATH_MSG_INFO( m_layerRecords.size() << " LayerMaterialRecords to be finalized for this material mapping run." );
-
+ ATH_MSG_INFO("========================================================================================= ");
+ ATH_MSG_INFO("finalize() starts ...");
+
+#ifdef TRKDETDESCR_MEMUSAGE
+ m_memoryLogger.refresh(getpid());
+ ATH_MSG_INFO("[ memory usage ] Start building of material maps: " );
+ ATH_MSG_INFO( m_memoryLogger );
+#endif
+
+ // create a dedicated LayerMaterialMap by layerMaterialCreator;
+ std::map< std::string, Trk::LayerMaterialMap* > layerMaterialMaps;
+ for ( auto& lmcIter : m_layerMaterialCreators ){
+ ATH_MSG_INFO("-> Creating material map '"<< lmcIter->layerMaterialName() << "' from creator "<< lmcIter.typeAndName() );
+ layerMaterialMaps[lmcIter->layerMaterialName()] = new Trk::LayerMaterialMap();
+ }
+
+ ATH_MSG_INFO( m_layerRecords.size() << "LayerMaterialRecords to be finalized for this material mapping run.");
+
// loop over the layers and output the stuff --- fill the associatedLayerMaterialProperties
- auto lIter = m_layerRecords.begin();
- for ( ; lIter != m_layerRecords.end(); ++lIter) {
+ for ( auto& lIter : m_layerRecords ) {
// Get the key map, the layer & the volume name
- const Trk::Layer* layer = (lIter->first);
- Trk::LayerIndex layerKey = (layer->layerIndex());
+ const Trk::Layer* layer = lIter.first;
+ Trk::LayerIndex layerKey = layer->layerIndex();
// get the enclosing tracking volume
const Trk::TrackingVolume* eVolume = layer->enclosingTrackingVolume();
// assign the string
- std::string vName = eVolume ? (eVolume->volumeName()) : "Unknown";
-
- ATH_MSG_INFO( "Finalize MaterialAssociation for Layer " << layerKey.value() << " in " << eVolume->volumeName() );
+ std::string vName = eVolume ? (eVolume->volumeName()) : " BoundaryCollection ";
+ ATH_MSG_INFO("Finalize MaterialAssociation for Layer "<< layerKey.value() << " in " << vName );
// normalize - use m_finalizeRunDebug
- (lIter->second).finalizeRun(m_layerMaterialScreenOutput);
-
- // the new ones to create
- Trk::LayerMaterialProperties* layerMaterialProperties = 0;
- // get the vector of MaterialProperties
- const Trk::MaterialPropertiesMatrix& recordedLayerMaterialProperties = (lIter->second).associatedLayerMaterial();
- if (recordedLayerMaterialProperties.size()==1 && recordedLayerMaterialProperties[0].size()==1) {
- // create homogeneous
- const Trk::MaterialProperties* singleProperties = (recordedLayerMaterialProperties[0][0]) ?
- (recordedLayerMaterialProperties[0][0]) : new Trk::MaterialProperties();
- layerMaterialProperties = new Trk::HomogeneousLayerMaterial(*singleProperties,1.);
- delete singleProperties;
- } else if (recordedLayerMaterialProperties.size()) {
- // get the BinUtility
- const Trk::BinUtility& recordLayerMaterialBinUtil = (*(lIter->second).binUtility());
- // create the new BinnedMaterial
- Trk::BinnedLayerMaterial* binnedMaterial = new Trk::BinnedLayerMaterial(recordLayerMaterialBinUtil,
- recordedLayerMaterialProperties,
- 1.);
- // the compress option for 2D material - simple gif algorithm
- if (m_compressMaterial && recordLayerMaterialBinUtil.max(1) > 1) {
- layerMaterialProperties = compressMaterial(*binnedMaterial);
- delete binnedMaterial;
- } else // this is the option for 1D -> don't compress
- layerMaterialProperties = binnedMaterial;
+ (lIter.second).finalizeRun(m_mapComposition);
+ // output the material to the analyser if registered
+ if (!m_layerMaterialRecordAnalyser.empty() && m_layerMaterialRecordAnalyser->analyseLayerMaterial(*layer, lIter.second).isFailure() )
+ ATH_MSG_WARNING("Could not analyse the LayerMaterialRecord for layer "<< layerKey.value() );
+ // check if we have analysers per creator
+ bool analyse = (m_layerMaterialCreators.size() == m_layerMaterialAnalysers.size());
+ // and now use the creators to make the maps out of the LayerMaterialRecord
+ size_t ilmc = 0;
+ for ( auto& lmcIter : m_layerMaterialCreators ){
+ // call the creator and register in the according map
+#ifdef TRKDETDESCR_MEMUSAGE
+ m_memoryLogger.refresh(getpid());
+ ATH_MSG_INFO("[ memory usage ] Before building the map for Layer "<< layerKey.value() );
+ ATH_MSG_INFO( m_memoryLogger );
+#endif
+ const Trk::LayerMaterialProperties* lMaterial = lmcIter->createLayerMaterial(*layer, lIter.second);
+#ifdef TRKDETDESCR_MEMUSAGE
+ m_memoryLogger.refresh(getpid());
+ ATH_MSG_INFO("[ memory usage ] After building the map for Layer "<< layerKey.value() );
+ ATH_MSG_INFO( m_memoryLogger );
+#endif
+ if (lMaterial)
+ ATH_MSG_VERBOSE("LayerMaterial map created as "<< *lMaterial );
+ // insert the created map for the given layer
+ (*layerMaterialMaps[lmcIter->layerMaterialName()])[layerKey.value()] = lMaterial;
+ // analyse the it if configured
+ if (analyse && lMaterial && (m_layerMaterialAnalysers[ilmc]->analyseLayerMaterial(*layer, *lMaterial)).isFailure() )
+ ATH_MSG_WARNING("Could not analyse created LayerMaterialProperties for layer "<< layerKey.value() );
+ ++ilmc;
}
- // fill them into the map that is written to POOL
- (*associatedLayerMaterialProperties)[layerKey] = layerMaterialProperties;
}
- ATH_MSG_INFO( " LayerMaterialProperties for " << associatedLayerMaterialProperties->size() << " layers recorded." );
-
- // write output to detector store
- if ( (detStore()->record(associatedLayerMaterialProperties, m_outputLayerMaterialSetName, false)).isFailure())
- ATH_MSG_ERROR( "Writing to DetectorStore was not successful!" );
- else {
- ATH_MSG_INFO( "LayerMaterialPropertiesMap: " << m_outputLayerMaterialSetName << " written to the DetectorStore!" );
+ ATH_MSG_INFO("Finalize map synchronization and write the maps to the DetectorStore.");
+
+ for (auto& ilmIter : layerMaterialMaps ){
+ // elementTable handling - if existent
+ if (m_mapComposition){
+ Trk::SharedObject<const Trk::ElementTable> tElementTable(new Trk::ElementTable(*m_elementTable));
+ ilmIter.second->updateElementTable(tElementTable);
+ if (ilmIter.second->elementTable()){
+ ATH_MSG_INFO("ElementTable for LayerMaterialMap '" << ilmIter.first << "' found and syncrhonized." );
+ ATH_MSG_INFO( *(ilmIter.second->elementTable()) );
+ }
+ }
+ // detector store writing
+ if ( (detStore()->record(ilmIter.second, ilmIter.first, false)).isFailure()){
+ ATH_MSG_ERROR( "Writing of LayerMaterialMap with name '" << ilmIter.first << "' was not successful." );
+ delete ilmIter.second;
+ } else ATH_MSG_INFO( "LayerMaterialMap: " << ilmIter.first << " written to the DetectorStore!" );
}
+ delete m_elementTable;
+
+#ifdef TRKDETDESCR_MEMUSAGE
+ m_memoryLogger.refresh(getpid());
+ ATH_MSG_INFO( "[ memory usage ] At the end of the material map creation.");
+ ATH_MSG_INFO( m_memoryLogger );
+#endif
ATH_MSG_INFO( "========================================================================================= " );
ATH_MSG_INFO( " -> Total mapped hits : " << m_mapped );
double unmapped = (m_unmapped+m_mapped) ? double(m_unmapped)/double(m_unmapped+m_mapped) : 0.;
ATH_MSG_INFO( " -> Total (rel.) unmapped hits : " << m_unmapped << " (" << unmapped << ")" );
- ATH_MSG_INFO(" -> Skipped (outisde) : " << m_skippedOutside );
+ ATH_MSG_INFO( " -> Skipped (outisde) : " << m_skippedOutside );
ATH_MSG_INFO( "========================================================================================= " );
ATH_MSG_INFO( "finalize() successful");
-
return StatusCode::SUCCESS;
-
}
@@ -454,19 +388,28 @@ StatusCode Trk::MaterialMapping::retrieveTrackingGeometry()
// Retrieve the TrackingGeometry from the DetectorStore
if ((detStore()->retrieve(m_trackingGeometry, m_trackingGeometrySvc->trackingGeometryName())).isFailure()) {
- ATH_MSG_FATAL( "Could not retrieve TrackingGeometry from DetectorStore!" );
+ ATH_MSG_FATAL("Could not retrieve TrackingGeometry from DetectorStore!");
return StatusCode::FAILURE;
}
+
// either get a string volume or the highest one
-
const Trk::TrackingVolume* trackingVolume = m_trackingGeometry->highestTrackingVolume();
// prepare the mapping volume
m_mappingVolume = m_trackingGeometry->trackingVolume(m_mappingVolumeName);
+ // register the confined layers from the TrackingVolume
registerVolume(*trackingVolume, 0);
+
+ ATH_MSG_INFO("Add "<< m_layerRecords.size() << " confined volume layers to mapping setup.");
+ ATH_MSG_INFO("Add "<< m_trackingGeometry->boundaryLayers().size() << " boundary layers to mapping setup.");
+
+ // register the layers from boundary surfaces
+ auto bLayerIter = m_trackingGeometry->boundaryLayers().begin();
+ for (; bLayerIter != m_trackingGeometry->boundaryLayers().end(); ++bLayerIter)
+ insertLayerMaterialRecord(*bLayerIter->first);
- ATH_MSG_INFO( "Map for " << m_layerRecords.size() << " Layers booked & prepared" );
+ ATH_MSG_INFO("Map for "<< m_layerRecords.size() << " layers booked & prepared for mapping procedure");
return StatusCode::SUCCESS;
@@ -478,8 +421,8 @@ void Trk::MaterialMapping::registerVolume(const Trk::TrackingVolume& tvol, int l
int sublevel = lvl+1;
for (int indent=0; indent<sublevel; ++indent)
- std::cout << " ";
- std::cout << "TrackingVolume name: " << tvol.volumeName() << std::endl;
+ std::cout << " ";
+ std::cout << "TrackingVolume name: "<< tvol.volumeName() << std::endl;
// all those to be processed
std::vector<const Trk::Layer*> volumeLayers;
@@ -490,8 +433,8 @@ void Trk::MaterialMapping::registerVolume(const Trk::TrackingVolume& tvol, int l
// this go ahead with the layers
const std::vector<const Trk::Layer*>& layers = confinedLayers->arrayObjects();
for (int indent=0; indent<sublevel; ++indent)
- std::cout << " ";
- std::cout << "- found : " << layers.size() << " confined Layers" << std::endl;
+ std::cout << " ";
+ std::cout << "- found : "<< layers.size() << "confined Layers"<< std::endl;
// loop over and fill them
auto clIter = layers.begin();
auto clIterE = layers.end();
@@ -503,58 +446,19 @@ void Trk::MaterialMapping::registerVolume(const Trk::TrackingVolume& tvol, int l
volumeLayers.push_back((*clIter));
}
}
-
- // collect all entry layers treatment
- const Trk::EntryLayerProvider* entryLayerProvider = tvol.entryLayerProvider();
- if (entryLayerProvider) {
- // -------------------------
- const std::vector<const Trk::Layer*>& entryLayers = entryLayerProvider->layers();
- for (int indent=0; indent<sublevel; ++indent)
- std::cout << " ";
- std::cout << "- found : " << entryLayers.size() << " entry Layers" << std::endl;
- auto elIter = entryLayers.begin();
- auto elIterE = entryLayers.end();
- for ( ; elIter != elIterE; ++elIter ) {
- const Amg::Vector3D& sReferencePoint = (*elIter)->surfaceRepresentation().globalReferencePoint();
- bool insideMappingVolume = m_mappingVolume ? m_mappingVolume->inside(sReferencePoint) : true;
- if ((*elIter)->layerMaterialProperties() && insideMappingVolume)
- volumeLayers.push_back((*elIter));
- }
- }
-
// now create LayerMaterialRecords for all
- auto lIter = volumeLayers.begin();
- auto lIterE = volumeLayers.end();
- for ( ; lIter != lIterE; ++lIter ){
- // first occurrance, create a new LayerMaterialRecord
- // - get the bin utility for the binned material (if necessary)
- // - get the material first
- const Trk::LayerMaterialProperties* layerMaterialProperties = (*lIter)->layerMaterialProperties();
- // - dynamic cast to the BinnedLayerMaterial
- const Trk::BinnedLayerMaterial* layerBinnedMaterial
- = dynamic_cast<const Trk::BinnedLayerMaterial*>(layerMaterialProperties);
- // get the binned array
- const Trk::BinUtility* layerMaterialBinUtility = (layerBinnedMaterial) ? layerBinnedMaterial->binUtility() : 0;
- // now fill the layer material record
- if (layerMaterialBinUtility){
- // create a new Layer Material record in the map
- Trk::LayerMaterialRecord lmr((*lIter)->thickness(),
- layerMaterialBinUtility,
- (Trk::MaterialAssociationType)m_associationType);
- // and fill it into the map
- m_layerRecords[(*lIter)] = lmr;
- }
- }
-
+ for ( auto& lIter : volumeLayers )
+ insertLayerMaterialRecord(*lIter);
+
// step dopwn the navigation tree to reach the confined volumes
const Trk::BinnedArray< Trk::TrackingVolume >* confinedVolumes = tvol.confinedVolumes();
if (confinedVolumes) {
const std::vector<const Trk::TrackingVolume*>& volumes = confinedVolumes->arrayObjects();
for (int indent=0; indent<sublevel; ++indent)
- std::cout << " ";
- std::cout << "- found : " << volumes.size() << " confined TrackingVolumes" << std::endl;
+ std::cout << " ";
+ std::cout << "- found : "<< volumes.size() << "confined TrackingVolumes"<< std::endl;
// loop over the confined volumes
auto volumesIter = volumes.begin();
for (; volumesIter != volumes.end(); ++volumesIter)
@@ -565,146 +469,25 @@ void Trk::MaterialMapping::registerVolume(const Trk::TrackingVolume& tvol, int l
}
-Trk::CompressedLayerMaterial* Trk::MaterialMapping::compressMaterial(const Trk::BinnedLayerMaterial& binnedMaterial) {
- // get the matrix
- const Trk::MaterialPropertiesMatrix& materialProperties = binnedMaterial.fullMaterial();
- // the vector to be created and reserve the maximum
- Trk::MaterialPropertiesVector materialVector;
- materialVector.reserve(m_compressedMaterialZARhoBins*m_compressedMaterialX0Bins);
- // this method creates a compressed representation of the BinnedLayerMaterial
- const Trk::BinUtility* cbinutil = binnedMaterial.binUtility();
- // nF x nS
- size_t nFirstBins = cbinutil->max(0)+1;
- size_t nSecondBins = cbinutil->max(1)+1;
- // low, high boundaries
- double x0min = 10e10;
- double x0max = 0.;
- double avZArhoMin = 10e10;
- double avZArhoMax = 0.;
- // create two maps, the compression map and the index map
- std::vector< std::vector<unsigned short int> > materialBins;
- // (1) FIRST LOOP, get boundaries
- materialBins.reserve(nSecondBins);
- for (size_t isec = 0; isec < nSecondBins; ++isec) {
- std::vector<unsigned short int> firstbins(nFirstBins,0);
- materialBins.push_back(firstbins);
- // loop over the bins
- for (size_t ifir = 0; ifir < nFirstBins; ++ifir) {
- // get the current material properties
- const Trk::MaterialProperties* matProp = materialProperties[isec][ifir];
- if (matProp) {
- double tinX0 = matProp->thicknessInX0();
- double avZArho = matProp->zOverAtimesRho();
- x0min = tinX0 < x0min ? tinX0 : x0min;
- x0max = tinX0 > x0max ? tinX0 : x0max;
- avZArhoMin = avZArho < avZArhoMin ? avZArho : avZArhoMin;
- avZArhoMax = avZArho > avZArhoMax ? avZArho : avZArhoMax;
- }
- }
- }
- // min / max is defined, find step size
- double stepX0 = (x0max-x0min)/m_compressedMaterialX0Bins;
- double stepZArho = (avZArhoMax-avZArhoMin)/m_compressedMaterialZARhoBins;
- // get the material histogram
- std::vector< std::vector< std::vector< Trk::IndexedMaterial> > > materialHistogram;
- materialHistogram.reserve(m_compressedMaterialZARhoBins);
- // prepare the histogram
- for (size_t izarho = 0; izarho < m_compressedMaterialZARhoBins; ++izarho) {
- std::vector< std::vector < Trk::IndexedMaterial > > x0materialbins;
- x0materialbins.reserve(m_compressedMaterialX0Bins);
- for (size_t ix0 = 0; ix0 < m_compressedMaterialX0Bins; ++ix0) {
- std::vector < Trk::IndexedMaterial > materialBin;
- x0materialbins.push_back( materialBin );
- }
- materialHistogram.push_back(x0materialbins);
- }
- // fill the histogram
- for (size_t isec = 0; isec < nSecondBins; ++isec) {
- for (size_t ifir = 0; ifir < nFirstBins; ++ifir) {
- // get the material properties
- const Trk::MaterialProperties* matProp = dynamic_cast<const Trk::MaterialProperties*>(materialProperties[isec][ifir]);
- if (matProp) {
- // calculate the bins of the material histogram
- double tinX0 = matProp->thicknessInX0();
- double avZArho = matProp->zOverAtimesRho();
- int x0bin = int( (tinX0-x0min)/stepX0 );
- int zarhobin = int( (avZArho-avZArhoMin)/stepZArho );
- // range protection
- x0bin = ( (size_t)x0bin >= m_compressedMaterialX0Bins) ? m_compressedMaterialX0Bins-1 : x0bin;
- x0bin = x0bin < 0 ? 0 : x0bin;
- zarhobin = ( (size_t)zarhobin >= m_compressedMaterialZARhoBins) ? m_compressedMaterialZARhoBins-1 : zarhobin;
- zarhobin = zarhobin < 0 ? 0 : zarhobin;
- // create indexed material
- Trk::IndexedMaterial idxMaterial;
- idxMaterial.materialProperties = matProp;
- idxMaterial.firstBin = ifir;
- idxMaterial.secondBin = isec;
- // fill into the material histogram
- materialHistogram[zarhobin][x0bin].push_back(idxMaterial);
- }
- }
- }
- // merge the bins and ready
- materialVector.push_back(0);
- // prepare the histogram
- for (size_t izarho = 0; izarho < m_compressedMaterialZARhoBins; ++izarho) {
- for (size_t ix0 = 0; ix0 < m_compressedMaterialX0Bins; ++ix0) {
- // get the indexed material properties
- std::vector< Trk::IndexedMaterial > indexedMaterial = materialHistogram[izarho][ix0];
- if (indexedMaterial.size()) {
- double avT = 0.; // thickness: by default on one layer it should be the same !
- double tinX0 = 0.;
- double tinL0 = 0.;
- double avA = 0.;
- double avZ = 0.;
- double avRho = 0.;
- std::vector< Trk::IndexedMaterial >::iterator idmIter = indexedMaterial.begin();
- std::vector< Trk::IndexedMaterial >::iterator idmIterEnd = indexedMaterial.end();
- for ( ; idmIter != idmIterEnd; ++idmIter ) {
- tinX0 += (*idmIter).materialProperties->thicknessInX0();
- tinL0 += (*idmIter).materialProperties->thicknessInL0();
- avA += (*idmIter).materialProperties->averageA();
- avZ += (*idmIter).materialProperties->averageZ();
- avRho += (*idmIter).materialProperties->averageRho();
- }
- double measure = 1./(indexedMaterial.size());
- // average it
- tinX0 *= measure;
- tinL0 *= measure;
- avA *= measure;
- avZ *= measure;
- avRho *= measure;
- avT *= measure;
- // compress to a model thickness [ rho affected ]
- avRho *= avT/m_compressedMaterialThickness;
- materialVector.push_back(new Trk::MaterialProperties(m_compressedMaterialThickness,
- m_compressedMaterialThickness/tinX0,
- m_compressedMaterialThickness/tinL0,
- avA,
- avZ,
- avRho));
- // now set the index
- int matindex = int(materialVector.size()-1);
- idmIter = indexedMaterial.begin();
- for ( ; idmIter != idmIterEnd; ++idmIter )
- materialBins[(*idmIter).secondBin][(*idmIter).firstBin] = matindex;
- }
- }
- }
-
- // change the 2bin matrix to a 1bin vector (better for persistency)
- std::vector<unsigned short int> materialBinsVector;
- materialBinsVector.reserve( (cbinutil->max(0)+1)*(cbinutil->max(1)+1) );
- std::vector< std::vector<unsigned short int> >::iterator binVecIter = materialBins.begin();
- std::vector< std::vector<unsigned short int> >::iterator binVecIterEnd = materialBins.end();
- for ( ; binVecIter != binVecIterEnd; ++binVecIter) {
- std::vector<unsigned short int>::iterator binIter = (*binVecIter).begin();
- std::vector<unsigned short int>::iterator binIterEnd = (*binVecIter).end();
- for ( ; binIter != binIterEnd; ++binIter )
- materialBinsVector.push_back(*binIter);
- }
- // create the compressed material
- return new Trk::CompressedLayerMaterial(*cbinutil,materialVector,materialBinsVector);
+void Trk::MaterialMapping::insertLayerMaterialRecord(const Trk::Layer& lay){
+ // first occurrance, create a new LayerMaterialRecord
+ // - get the bin utility for the binned material (if necessary)
+ // - get the material first
+ const Trk::LayerMaterialProperties* layerMaterialProperties = lay.layerMaterialProperties();
+ // - dynamic cast to the BinnedLayerMaterial
+ const Trk::BinnedLayerMaterial* layerBinnedMaterial
+ = dynamic_cast<const Trk::BinnedLayerMaterial*>(layerMaterialProperties);
+ // get the binned array
+ const Trk::BinUtility* layerMaterialBinUtility = (layerBinnedMaterial) ? layerBinnedMaterial->binUtility() : 0;
+ // now fill the layer material record
+ if (layerMaterialBinUtility){
+ // create a new Layer Material record in the map
+ Trk::LayerMaterialRecord lmr((m_useLayerThickness ? lay.thickness() : 1.),
+ layerMaterialBinUtility,
+ (Trk::MaterialAssociationType)m_associationType);
+ // and fill it into the map
+ m_layerRecords[&lay] = lmr;
+ }
}
diff --git a/Tracking/TrkDetDescr/TrkDetDescrAlgs/src/MaterialValidation.cxx b/Tracking/TrkDetDescr/TrkDetDescrAlgs/src/MaterialValidation.cxx
index 432e77762234b193a69f363ad0ba7d6e91131fa2..102d5e9e105bf3ba1b230a8cd5d8500039d9325e 100755
--- a/Tracking/TrkDetDescr/TrkDetDescrAlgs/src/MaterialValidation.cxx
+++ b/Tracking/TrkDetDescr/TrkDetDescrAlgs/src/MaterialValidation.cxx
@@ -10,7 +10,6 @@
#include "GeoPrimitives/GeoPrimitivesToStringConverter.h"
// Trk
#include "TrkDetDescrAlgs/MaterialValidation.h"
-#include "TrkGeometry/EntryLayerProvider.h"
#include "TrkGeometry/TrackingVolume.h"
#include "TrkGeometry/TrackingGeometry.h"
#include "TrkGeometry/Layer.h"
@@ -19,6 +18,7 @@
#include "TrkGeometry/AssociatedMaterial.h"
#include "TrkDetDescrInterfaces/ITrackingGeometrySvc.h"
#include "TrkDetDescrInterfaces/IMaterialMapper.h"
+#include "TrkNeutralParameters/NeutralParameters.h"
// test
#include "TrkVolumes/CylinderVolumeBounds.h"
@@ -36,7 +36,9 @@ Trk::MaterialValidation::MaterialValidation(const std::string& name, ISvcLocator
m_maxMaterialMappingEvents(25000),
m_flatDist(0),
m_etaMin(-3.),
- m_etaMax(3.)
+ m_etaMax(3.),
+ m_runNativeNavigation(true),
+ m_accTinX0(0)
{
// ---------------------- The TrackingGeometrySvc ------------------------ //
@@ -48,6 +50,8 @@ Trk::MaterialValidation::MaterialValidation(const std::string& name, ISvcLocator
// ---------------------- Range setup ----------------------------------- //
declareProperty("MinEta" , m_etaMin);
declareProperty("MaxEta" , m_etaMax);
+ // ---------------------- Native navigation ----------------------------- //
+ declareProperty("NativeNavigation" , m_runNativeNavigation);
}
Trk::MaterialValidation::~MaterialValidation()
@@ -77,7 +81,7 @@ StatusCode Trk::MaterialValidation::initialize()
StatusCode Trk::MaterialValidation::execute()
{
- ATH_MSG_VERBOSE( "MaterialValidation execute() start" );
+ ATH_MSG_VERBOSE( "MaterialValidation execute() start ================================================" );
// ------------------------------- get the trackingGeometry at first place
if (!m_trackingGeometry) {
@@ -92,10 +96,13 @@ StatusCode Trk::MaterialValidation::execute()
double eta = m_etaMin + (m_etaMax-m_etaMin)*m_flatDist->shoot();
double theta = 2.*atan(exp(-eta));
double phi = M_PI * ( 2*m_flatDist->shoot() - 1.);
+ m_accTinX0 = 0;
// get the position and riection from the random numbers
Amg::Vector3D position(0.,0.,0.);
Amg::Vector3D direction(sin(theta)*cos(phi),sin(theta)*sin(phi), cos(theta));
+
+ ATH_MSG_DEBUG("[>] Start mapping event with phi | eta = " << phi << " | " << direction.eta());
// find the start TrackingVolume
const Trk::TrackingVolume* sVolume = m_trackingGeometry->lowestTrackingVolume(position);
@@ -104,11 +111,8 @@ StatusCode Trk::MaterialValidation::execute()
Trk::PositionAtBoundary paB = collectMaterialAndExit(*nVolume, position, direction);
position = paB.first;
nVolume = paB.second;
- }
-
- // finalize the event
- m_materialMapper->finalizeEvent(eta,phi);
-
+ }
+ ATH_MSG_DEBUG("[<] Finishing event with collected path [X0] = " << m_accTinX0);
return StatusCode::SUCCESS;
}
@@ -117,135 +121,172 @@ Trk::PositionAtBoundary Trk::MaterialValidation::collectMaterialAndExit(const Tr
const Amg::Vector3D& direction)
{
// get the entry layers -----------------------------------------------------------
- std::map<double, std::pair<const Trk::Layer*, Amg::Vector3D> > intersectedLayers;
std::map<double, Trk::AssociatedMaterial> collectedMaterial;
- ATH_MSG_VERBOSE("Entering Volume: " << tvol.volumeName() << "- at " << Amg::toString(position));
+ // all boundaries found --- proceed
+ Trk::PositionAtBoundary pab(position, 0);
- // get theta from the direction
- double theta = direction.theta();
+ ATH_MSG_DEBUG("[>>] Entering Volume: " << tvol.volumeName() << "- at " << Amg::toString(position));
- // Process the entry layers if they exist
- const Trk::EntryLayerProvider* entryLayerProvider = tvol.entryLayerProvider();
- if (entryLayerProvider) {
- // display output
- const std::vector<const Trk::Layer*>& entryLayers = entryLayerProvider->layers();
- auto eLayIter = entryLayers.begin();
- auto eLayIterE = entryLayers.end();
- for ( ; eLayIter != eLayIterE; ++eLayIter){
- Trk::SurfaceIntersection esIntersection = (*eLayIter)->surfaceRepresentation().straightLineIntersection(position, direction, true, true);
- if (esIntersection.valid){
- // intersection valid : record it
- intersectedLayers[esIntersection.pathLength] = std::pair<const Trk::Layer*, Amg::Vector3D>(*eLayIter,esIntersection.intersection);
- // calculate the pathCorrection
- double stepLength = (*eLayIter)->surfaceRepresentation().type() == Trk::Surface::Cylinder ? fabs(1./sin(theta)) : fabs(1./cos(theta));
- // get & record the material
- const Trk::MaterialProperties* mprop = 0;
- if ((*eLayIter)->layerMaterialProperties()){
- Amg::Vector3D mposition = (*eLayIter)->surfaceRepresentation().transform().inverse()*esIntersection.intersection;
- mprop = (*eLayIter)->layerMaterialProperties()->fullMaterial(mposition);
- }
- // if you have material properties
- if (mprop)
- collectedMaterial[esIntersection.pathLength] = Trk::AssociatedMaterial(esIntersection.intersection, mprop, stepLength, &tvol, (*eLayIter));
- else
- collectedMaterial[esIntersection.pathLength] = Trk::AssociatedMaterial(esIntersection.intersection, &tvol, (*eLayIter));
-
- ATH_MSG_VERBOSE(" - record material hit at (entry) layer with index " << (*eLayIter)->layerIndex().value() << " - at " << Amg::toString(esIntersection.intersection) );
- if (mprop)
- ATH_MSG_VERBOSE(" - MaterialProperties are " << (*mprop) );
- else
- ATH_MSG_VERBOSE(" - No MaterialProperties found." );
-
- }
- }
- }
+ Trk::NeutralCurvilinearParameters cvp(position,direction,0.);
- // Process the contained layers if they exist
- const Trk::LayerArray* layerArray = tvol.confinedLayers();
- if (layerArray) {
- // display output
- const std::vector<const Trk::Layer*>& layers = layerArray->arrayObjects();
- auto layIter = layers.begin();
- auto layIterE = layers.end();
- for ( ; layIter != layIterE; ++layIter){
- if ( (*layIter)->layerMaterialProperties() ){
- Trk::SurfaceIntersection lsIntersection = (*layIter)->surfaceRepresentation().straightLineIntersection(position, direction, true, true);
- if (lsIntersection.valid){
- intersectedLayers[lsIntersection.pathLength] = std::pair<const Trk::Layer*, Amg::Vector3D>(*layIter,lsIntersection.intersection);
- // calculate the pathCorrection
- double stepLength = (*layIter)->surfaceRepresentation().type() == Trk::Surface::Cylinder ? fabs(1./sin(theta)) : fabs(1./cos(theta));
- // get & record the material
- // - the position on the surface
- Amg::Vector3D mposition = (*layIter)->surfaceRepresentation().transform().inverse()*lsIntersection.intersection;
- const Trk::MaterialProperties* mprop = (*layIter)->layerMaterialProperties()->fullMaterial(mposition);
- if (mprop)
- collectedMaterial[lsIntersection.pathLength] = Trk::AssociatedMaterial(lsIntersection.intersection, mprop, stepLength, &tvol, (*layIter));
- else
- collectedMaterial[lsIntersection.pathLength] = Trk::AssociatedMaterial(lsIntersection.intersection, &tvol, (*layIter));
-
- ATH_MSG_VERBOSE(" - record material hit at (conf.) layer with index " << (*layIter)->layerIndex().value() << " - at " << Amg::toString(lsIntersection.intersection) );
- if (mprop)
- ATH_MSG_VERBOSE(" - MaterialProperties are " << (*mprop) );
- else
- ATH_MSG_VERBOSE(" - No MaterialProperties found." );
+ if (m_runNativeNavigation){
+ // A : collect all hit layers
+ auto layerIntersections = tvol.materialLayersOrdered<Trk::NeutralCurvilinearParameters>(NULL,NULL,cvp,Trk::alongMomentum);
+ // loop over the layers
+ for (auto& lCandidate : layerIntersections ) {
+ // get the layer
+ const Trk::Layer* layer = lCandidate.object;
+ double pathLength = lCandidate.intersection.pathLength;
+ // get the associate material
+ if (layer->layerMaterialProperties()){
+ // take it from the global position
+ const Trk::MaterialProperties* mprop = layer->layerMaterialProperties()->fullMaterial(lCandidate.intersection.position);
+ if (mprop){
+ double stepLength = mprop->thickness()*fabs(layer->surfaceRepresentation().pathCorrection(lCandidate.intersection.position,direction));
+ collectedMaterial[pathLength] = Trk::AssociatedMaterial(lCandidate.intersection.position, mprop, stepLength, &tvol, layer);
+ }
+ }
+ }
+ // B : collect all boundary layers, start from last hit layer
+ Amg::Vector3D lastPosition = collectedMaterial.size() ? collectedMaterial.rbegin()->second.materialPosition() : (position + direction.unit());
+ Trk::NeutralCurvilinearParameters lcp(lastPosition,direction,0.);
+ // boundary surfaces
+ auto boundaryIntersections = tvol.boundarySurfacesOrdered<Trk::NeutralCurvilinearParameters>(lcp,Trk::alongMomentum);
+ if (boundaryIntersections.size()){
+ // by definition is the first one
+ lastPosition = boundaryIntersections.begin()->intersection.position;
+ const Trk::BoundarySurface<Trk::TrackingVolume>* bSurfaceTV = boundaryIntersections.begin()->object;
+ const Trk::Surface& bSurface = bSurfaceTV->surfaceRepresentation();
+ // get the path lenght to it
+ if (bSurface.materialLayer() && bSurface.materialLayer()->layerMaterialProperties()){
+ const Trk::MaterialProperties* mprop = bSurface.materialLayer()->layerMaterialProperties()->fullMaterial(lastPosition);
+ double pathLength = (lastPosition-position).mag();
+ if (mprop){
+ double stepLength = mprop->thickness()*fabs(bSurface.pathCorrection(lastPosition,direction));
+ collectedMaterial[pathLength] = Trk::AssociatedMaterial(lastPosition, mprop, stepLength, &tvol, bSurface.materialLayer());
+ } else
+ collectedMaterial[pathLength] = Trk::AssociatedMaterial(lastPosition, &tvol, bSurface.materialLayer());
+ }
+ // set the new volume
+ const Trk::TrackingVolume* naVolume = bSurfaceTV->attachedVolume(lastPosition, direction, Trk::alongMomentum);
+ pab = Trk::PositionAtBoundary(lastPosition,naVolume);
+ } else
+ pab = Trk::PositionAtBoundary(lastPosition,0);
+ } else {
+ std::map<double, std::pair<const Trk::Layer*, Amg::Vector3D> > intersectedLayers;
+
+ // Process the contained layers if they exist
+ const Trk::LayerArray* layerArray = tvol.confinedLayers();
+ if (layerArray) {
+ // display output
+ const std::vector<const Trk::Layer*>& layers = layerArray->arrayObjects();
+ auto layIter = layers.begin();
+ auto layIterE = layers.end();
+ for ( ; layIter != layIterE; ++layIter){
+ if ( (*layIter)->layerMaterialProperties() ){
+ Trk::Intersection lsIntersection = (*layIter)->surfaceRepresentation().straightLineIntersection(position, direction, true, true);
+ if (lsIntersection.valid){
+ intersectedLayers[lsIntersection.pathLength] = std::pair<const Trk::Layer*, Amg::Vector3D>(*layIter,lsIntersection.position);
+ // get & record the material
+ // - the position on the surface
+ Amg::Vector3D mposition = (*layIter)->surfaceRepresentation().transform().inverse()*lsIntersection.position;
+ const Trk::MaterialProperties* mprop = (*layIter)->layerMaterialProperties()->fullMaterial(mposition);
+ if (mprop) {
+ double stepLength = mprop->thickness()*fabs((*layIter)->surfaceRepresentation().pathCorrection(lsIntersection.position,direction));
+ collectedMaterial[lsIntersection.pathLength] = Trk::AssociatedMaterial(lsIntersection.position, mprop, stepLength, &tvol, (*layIter));
+ } else
+ collectedMaterial[lsIntersection.pathLength] = Trk::AssociatedMaterial(lsIntersection.position, &tvol, (*layIter));
+
+ ATH_MSG_VERBOSE("[>>>>] record material hit at layer with index " << (*layIter)->layerIndex().value() << " - at " << Amg::toString(lsIntersection.position) );
+ if (mprop)
+ ATH_MSG_VERBOSE("[>>>>] MaterialProperties are " << (*mprop) );
+ else
+ ATH_MSG_VERBOSE("[>>>>] No MaterialProperties found." );
+ }
}
}
- }
- }
+ }
+
+ // material for confined layers collected, now go to boundary
+
+ // update the position to the last one
+ Amg::Vector3D lastPosition = intersectedLayers.size() ? (*(--(intersectedLayers.end()))).second.second : position;
+
+ std::map<double, Trk::VolumeExit > volumeExits;
+ // now find the exit point
+ const std::vector< Trk::SharedObject<const Trk::BoundarySurface<Trk::TrackingVolume> > >& bSurfaces = tvol.boundarySurfaces();
+ auto bSurfIter = bSurfaces.begin();
+ auto bSurfIterE = bSurfaces.end();
+ for ( ; bSurfIter != bSurfIterE; ++bSurfIter){
+ // omit positions on the surface
+ if ( !(*bSurfIter)->surfaceRepresentation().isOnSurface(lastPosition, true, 0.1, 0.1) ){
+ Trk::Intersection evIntersection = (*bSurfIter)->surfaceRepresentation().straightLineIntersection(lastPosition, direction, true, true);
+ ATH_MSG_VERBOSE("[>>>>] boundary surface intersection / validity :" << Amg::toString(evIntersection.position) << " / " << evIntersection.valid);
+ ATH_MSG_VERBOSE("[>>>>] with path length = " << evIntersection.pathLength );
+ if (evIntersection.valid){
+ // next attached Tracking Volume
+ const Trk::TrackingVolume* naVolume = (*bSurfIter)->attachedVolume(evIntersection.position, direction, Trk::alongMomentum);
+ // put it into the map
+ volumeExits[evIntersection.pathLength] = Trk::VolumeExit(naVolume, (&(*bSurfIter)->surfaceRepresentation()), evIntersection.position);
+ // volume exit
+ ATH_MSG_VERBOSE("[>>>>] found volume exit - at " << Amg::toString(evIntersection.position) );
+ }
+ } else
+ ATH_MSG_VERBOSE("[>>>>] starting position is on surface ! " );
+ }
+ // prepare the boundary
+ if (volumeExits.size()){
+ // get the first entry in the map: closest next volume
+ VolumeExit closestVolumeExit = (*volumeExits.begin()).second;
+ // check if the volume exit boundary has material attached
+ const Trk::Surface* bSurface = closestVolumeExit.bSurface;
+ if ( bSurface && bSurface->materialLayer() && bSurface->materialLayer()->layerMaterialProperties()){
+ ATH_MSG_VERBOSE("[>>>>] The boundary surface has an associated layer, collect material from there");
+ const Trk::MaterialProperties* mprop = bSurface->materialLayer()->layerMaterialProperties()->fullMaterial(closestVolumeExit.vExit);
+ double pathToExit = (closestVolumeExit.vExit-lastPosition).mag();
+ if (mprop){
+ double stepLength = mprop->thickness()*fabs(bSurface->pathCorrection(closestVolumeExit.vExit,direction));
+ collectedMaterial[pathToExit] = Trk::AssociatedMaterial(closestVolumeExit.vExit, mprop, stepLength, &tvol, bSurface->materialLayer());
+ } else
+ collectedMaterial[pathToExit] = Trk::AssociatedMaterial(closestVolumeExit.vExit, &tvol, bSurface->materialLayer());
+ }
+ //
+ if (closestVolumeExit.nVolume != &tvol && closestVolumeExit.nVolume) {
+ ATH_MSG_VERBOSE("[>>>>] Next Volume: " << closestVolumeExit.nVolume->volumeName() << " - at " << Amg::toString(closestVolumeExit.vExit) );
+ // return for further navigation
+ pab = Trk::PositionAtBoundary(closestVolumeExit.vExit,closestVolumeExit.nVolume);
+ }
+ } else {
+ ATH_MSG_VERBOSE( "[>>>>] No exit found from Volume '" << tvol.volumeName() << "' - starting radius = " << lastPosition.perp() );
+ const Trk::CylinderVolumeBounds* cvb = dynamic_cast<const Trk::CylinderVolumeBounds*>(&(tvol.volumeBounds()));
+ if (cvb)
+ ATH_MSG_VERBOSE( "[>>>>] Volume outer radius = " << cvb->outerRadius() );
+ }
+ }
+ // finally collect the material
+ ATH_MSG_DEBUG("[>>>] Collecting materials from "<< collectedMaterial.size() << " layers");
// provide the material to the material mapper
auto cmIter = collectedMaterial.begin();
auto cmIterE = collectedMaterial.end();
- for ( ; cmIter != cmIterE; ++cmIter )
+ for ( ; cmIter != cmIterE; ++cmIter ){
m_materialMapper->recordMaterialHit(cmIter->second, cmIter->second.materialPosition());
-
- // update the position to the last one
- Amg::Vector3D lastPosition = intersectedLayers.size() ? (*(--(intersectedLayers.end()))).second.second : position;
- // exit the volume
- std::map<double, std::pair<const Trk::TrackingVolume*, Amg::Vector3D> > volumeExits;
- // now find the exit point
- const std::vector< Trk::SharedObject<const Trk::BoundarySurface<Trk::TrackingVolume> > >& bSurfaces = tvol.boundarySurfaces();
- auto bSurfIter = bSurfaces.begin();
- auto bSurfIterE = bSurfaces.end();
- for ( ; bSurfIter != bSurfIterE; ++bSurfIter){
- // omit positions on the surface
- if ( !(*bSurfIter)->surfaceRepresentation().isOnSurface(lastPosition, true, 0.1, 0.1) ){
- Trk::SurfaceIntersection evIntersection = (*bSurfIter)->surfaceRepresentation().straightLineIntersection(lastPosition, direction, true, true);
- ATH_MSG_DEBUG(" - boundary surface intersection / validity :" << Amg::toString(evIntersection.intersection) << " / " << evIntersection.valid);
- ATH_MSG_DEBUG(" - with path length = " << evIntersection.pathLength );
- if (evIntersection.valid){
- // next attached Tracking Volume
- const Trk::TrackingVolume* naVolume = (*bSurfIter)->attachedVolume(evIntersection.intersection, direction, Trk::alongMomentum);
- // put it into the map
- volumeExits[evIntersection.pathLength] = std::pair<const Trk::TrackingVolume*, Amg::Vector3D>(naVolume,evIntersection.intersection);
- // volume exit
- ATH_MSG_DEBUG(" - found volume exit - at " << Amg::toString(evIntersection.intersection) );
- }
- } else {
- ATH_MSG_DEBUG(" - starting position is on surface ! " );
- }
- }
-
- if (volumeExits.size()){
- // get the first entry in the map: closest next volume
- const Trk::TrackingVolume* cnVolume = (*volumeExits.begin()).second.first;
- if (cnVolume != &tvol) {
- if (cnVolume)
- ATH_MSG_VERBOSE("Next Volume: " << cnVolume->volumeName() << " - at " << Amg::toString((*volumeExits.begin()).second.second) );
- // record that
- m_materialMapper->record((*volumeExits.begin()).second.second);
- // return for further navigation
- return Trk::PositionAtBoundary((*volumeExits.begin()).second.second ,cnVolume);
+ m_accTinX0 += cmIter->second.steplengthInX0();
+ int layerIndex = cmIter->second.associatedLayer() ? cmIter->second.associatedLayer()->layerIndex().value() : 0;
+ ATH_MSG_DEBUG("[>>>] Accumulate pathLength/X0 on layer with index " << layerIndex << " - t/X0 (total so far) = " << cmIter->second.steplengthInX0() << " (" << m_accTinX0 << ")");
+ if (layerIndex){
+ std::string surfaceType = cmIter->second.associatedLayer()->surfaceRepresentation().type() == Trk::Surface::Cylinder ? "Cylinder at radius = " : "Disc at z-position = ";
+ std::string layerType = cmIter->second.associatedLayer()->surfaceArray() ? "Active " : "Passive ";
+ double rz = cmIter->second.associatedLayer()->surfaceRepresentation().type() == Trk::Surface::Cylinder ? cmIter->second.associatedLayer()->surfaceRepresentation().bounds().r() :
+ cmIter->second.associatedLayer()->surfaceRepresentation().center().z();
+ ATH_MSG_DEBUG(" " << layerType << surfaceType << rz);
}
- } else{
- ATH_MSG_DEBUG( "No exit found from Volume '" << tvol.volumeName() << "' - starting radius = " << lastPosition.perp() );
- const Trk::CylinderVolumeBounds* cvb = dynamic_cast<const Trk::CylinderVolumeBounds*>(&(tvol.volumeBounds()));
- if (cvb)
- ATH_MSG_DEBUG( " Volume outer radius = " << cvb->outerRadius() );
+ ATH_MSG_DEBUG(" Distance to origin is " << cmIter->second.materialPosition().mag() );
}
- return Trk::PositionAtBoundary(position, 0);
+ // return what you have
+ return pab;
}