DD4hep  1.32.1
Detector Description Toolkit for High Energy Physics
MaterialManager.cpp
Go to the documentation of this file.
1 //==========================================================================
2 // AIDA Detector description implementation
3 //--------------------------------------------------------------------------
4 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
5 // All rights reserved.
6 //
7 // For the licensing terms see $DD4hepINSTALL/LICENSE.
8 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
9 //
10 // Author : F.Gaede
11 //
12 //==========================================================================
13 #include "DDRec/MaterialManager.h"
14 #include "DD4hep/Exceptions.h"
15 #include "DD4hep/Detector.h"
16 
17 #include "TGeoVolume.h"
18 #include "TGeoManager.h"
19 #include "TGeoNode.h"
20 #include "TVirtualGeoTrack.h"
21 
22 #define MINSTEP 1.e-5
23 
24 namespace dd4hep {
25  namespace rec {
26 
27  MaterialManager::MaterialManager(Volume world) : _mV(0), _m( Material() ), _p0(),_p1(),_pos() {
28  _tgeoMgr = world->GetGeoManager();
29  }
30 
32 
33  }
34 
35  const PlacementVec& MaterialManager::placementsBetween(const Vector3D& p0, const Vector3D& p1 , double eps) {
36  materialsBetween(p0,p1,eps);
37  return _placeV;
38  }
39 
41  const Vector3D& p1,
42  double eps) {
43  const auto& materials = this->materialsBetween(p0, p1, eps);
44  return { materials, this->_placeV };
45  }
46 
47  const MaterialVec& MaterialManager::materialsBetween(const Vector3D& p0, const Vector3D& p1 , double eps) {
48  if( ( p0 != _p0 ) || ( p1 != _p1 ) ) {
49  // A backup is needed to restore the state of the navigator after the track is done
50  // see https://github.com/AIDASoft/DD4hep/issues/1413
51  _tgeoMgr->DoBackupState();
52  //---------------------------------------
53  _mV.clear() ;
54  _placeV.clear();
55  //
56  // algorithm copied from TGeoGearDistanceProperties.cc (A.Munnich):
57  //
58 
59  double startpoint[3], endpoint[3], direction[3];
60  double L=0;
61  for(unsigned int i=0; i<3; i++) {
62  startpoint[i] = p0[i];
63  endpoint[i] = p1[i];
64  direction[i] = endpoint[i] - startpoint[i];
65  L+=direction[i]*direction[i];
66  }
67  double totDist = sqrt( L ) ;
68 
69  //normalize direction
70  for(unsigned int i=0; i<3; i++)
71  direction[i]=direction[i]/totDist;
72 
73  _tgeoMgr->AddTrack(0, 12 ) ; // electron neutrino
74 
75  TGeoNode *node1 = _tgeoMgr->InitTrack(startpoint, direction);
76 
77  //check if there is a node at startpoint
78  if(!node1)
79  throw std::runtime_error("No geometry node found at given location. Either there is no node placed here or position is outside of top volume.");
80 
81  while ( !_tgeoMgr->IsOutside() ) {
82 
83  // TGeoNode *node2;
84  // TVirtualGeoTrack *track;
85 
86  // step to (and over) the next Boundary
87  TGeoNode * node2 = _tgeoMgr->FindNextBoundaryAndStep( 500, 1) ;
88 
89  if( !node2 || _tgeoMgr->IsOutside() )
90  break;
91 
92  const double *position = _tgeoMgr->GetCurrentPoint();
93  const double *previouspos = _tgeoMgr->GetLastPoint();
94 
95  double length = _tgeoMgr->GetStep();
96 
97  TVirtualGeoTrack *track = _tgeoMgr->GetLastTrack();
98 
99  //protection against infinitive loop in root which should not happen, but well it does...
100  //work around until solution within root can be found when the step gets very small e.g. 1e-10
101  //and the next boundary is never reached
102 
103 #if 1 //fg: is this still needed ?
104  if( length < MINSTEP ) {
105 
106  _tgeoMgr->SetCurrentPoint( position[0] + MINSTEP * direction[0],
107  position[1] + MINSTEP * direction[1],
108  position[2] + MINSTEP * direction[2] );
109 
110  length = _tgeoMgr->GetStep();
111  node2 = _tgeoMgr->FindNextBoundaryAndStep(500, 1) ;
112 
113  position = _tgeoMgr->GetCurrentPoint();
114  previouspos = _tgeoMgr->GetLastPoint();
115  }
116 #endif
117  // printf( " -- step length : %1.8e %1.8e %1.8e %1.8e %1.8e %1.8e %1.8e - %s \n" , length ,
118  // position[0], position[1], position[2], previouspos[0], previouspos[1], previouspos[2] , node1->GetMedium()->GetMaterial()->GetName() ) ;
119 
120  Vector3D posV( position ) ;
121 
122  double currDistance = ( posV - p0 ).r() ;
123 
124  // //if the next boundary is further than end point
125  // if(fabs(position[0])>fabs(endpoint[0]) || fabs(position[1])>fabs(endpoint[1])
126  // || fabs(position[2])>fabs(endpoint[2]))
127 
128  //if we travelled too far:
129  if( currDistance > totDist ) {
130 
131  length = sqrt( pow(endpoint[0]-previouspos[0],2) +
132  pow(endpoint[1]-previouspos[1],2) +
133  pow(endpoint[2]-previouspos[2],2) );
134 
135  track->AddPoint( endpoint[0], endpoint[1], endpoint[2], 0. );
136 
137 
138  if( length > eps ) {
139  _mV.emplace_back(node1->GetMedium(), length );
140  _placeV.emplace_back(node1,length);
141  }
142  break;
143  }
144 
145  track->AddPoint( position[0], position[1], position[2], 0.);
146 
147  if( length > eps ) {
148  _mV.emplace_back(node1->GetMedium(), length);
149  _placeV.emplace_back(node1,length);
150  }
151  node1 = node2;
152  }
153 
154 
155  //fg: protect against empty list:
156  if( _mV.empty() ){
157  _mV.emplace_back(node1->GetMedium(), totDist);
158  _placeV.emplace_back(node1,totDist);
159  }
160 
161 
162  _tgeoMgr->ClearTracks();
163 
164  _tgeoMgr->CleanGarbage();
165 
166  _tgeoMgr->DoRestoreState();
167 
168  //---------------------------------------
169 
170  _p0 = p0 ;
171  _p1 = p1 ;
172  }
173 
174  return _mV ;
175  }
176 
177 
179  if( pos != _pos ) {
180  TGeoNode *node = _tgeoMgr->FindNode( pos[0], pos[1], pos[2] ) ;
181  if( ! node ) {
182  std::stringstream err ;
183  err << " MaterialManager::material: No geometry node found at location: " << pos ;
184  throw std::runtime_error( err.str() );
185  }
186  _m = Material( node->GetMedium() );
187  _pv = node;
188  _pos = pos ;
189  }
190  return _m ;
191  }
192 
194  if( pos != _pos ) {
195  TGeoNode *node = _tgeoMgr->FindNode( pos[0], pos[1], pos[2] ) ;
196  if( ! node ) {
197  std::stringstream err ;
198  err << " MaterialManager::material: No geometry node found at location: " << pos ;
199  throw std::runtime_error( err.str() );
200  }
201  _m = Material( node->GetMedium() );
202  _pv = node;
203  _pos = pos;
204  }
205  return _pv;
206  }
207 
209 
210  std::stringstream sstr ;
211 
212  double sum_l = 0 ;
213  double sum_rho_l = 0 ;
214  double sum_rho_l_over_A = 0 ;
215  double sum_rho_l_Z_over_A = 0 ;
216  //double sum_rho_l_over_x = 0 ;
217  double sum_l_over_x = 0 ;
218  //double sum_rho_l_over_lambda = 0 ;
219  double sum_l_over_lambda = 0 ;
220 
221  for(unsigned i=0,n=materials.size(); i<n ; ++i){
222 
223  Material mat = materials[i].first ;
224  double l = materials[i].second ;
225 
226  if( i != 0 ) sstr << "_" ;
227  sstr << mat.name() << "_" << l ;
228 
229  double rho = mat.density() ;
230  double A = mat.A() ;
231  double Z = mat.Z() ;
232  double x = mat.radLength() ;
233  double lambda = mat.intLength() ;
234 
235  sum_l += l ;
236  sum_rho_l += rho * l ;
237  sum_rho_l_over_A += rho * l / A ;
238  sum_rho_l_Z_over_A += rho * l * Z / A ;
239  sum_l_over_x += l / x ;
240  sum_l_over_lambda += l / lambda ;
241  // sum_rho_l_over_x += rho * l / x ;
242  // sum_rho_l_over_lambda += rho * l / lambda ;
243  }
244 
245  double rho = sum_rho_l / sum_l ;
246 
247  double A = sum_rho_l / sum_rho_l_over_A ;
248  double Z = sum_rho_l_Z_over_A / sum_rho_l_over_A ;
249 
250  // radiation and interaction lengths already given in cm - average by length
251 
252  // double x = sum_rho_l / sum_rho_l_over_x ;
253  double x = sum_l / sum_l_over_x ;
254 
255  // double lambda = sum_rho_l / sum_rho_l_over_lambda ;
256  double lambda = sum_l / sum_l_over_lambda ;
257 
258 
259  return MaterialData( sstr.str() , Z, A, rho, x, lambda ) ;
260 
261  }
262 
263  } /* namespace rec */
264 } /* namespace dd4hep */
dd4hep::rec::MaterialManager::ScanData
Definition: MaterialManager.h:59
dd4hep::rec::MaterialManager::materialAt
const Material & materialAt(const Vector3D &pos)
Definition: MaterialManager.cpp:178
dd4hep::rec::MaterialManager::_mV
MaterialVec _mV
Cached materials.
Definition: MaterialManager.h:98
dd4hep::rec::MaterialManager::placementsBetween
const PlacementVec & placementsBetween(const Vector3D &p0, const Vector3D &p1, double eps=MaterialManager::epsilon)
Definition: MaterialManager.cpp:35
dd4hep::rec::Vector3D
Definition: Vector3D.h:32
dd4hep::rec::MaterialManager::materialsBetween
const MaterialVec & materialsBetween(const Vector3D &p0, const Vector3D &p1, double eps=MaterialManager::epsilon)
Definition: MaterialManager.cpp:47
Detector.h
dd4hep::rec::MaterialManager::createAveragedMaterial
MaterialData createAveragedMaterial(const MaterialVec &materials)
Definition: MaterialManager.cpp:208
dd4hep::PlacedVolume
Handle class holding a placed volume (also called physical volume)
Definition: Volumes.h:164
dd4hep::rec::MaterialManager::_p1
Vector3D _p1
Definition: MaterialManager.h:104
dd4hep::rec::MaterialVec
std::vector< std::pair< Material, double > > MaterialVec
Definition: MaterialManager.h:30
dd4hep::Handle::name
const char * name() const
Access the object name (or "" if not supported by the object)
dd4hep::Material::intLength
double intLength() const
Access the interaction length of the underlying material.
Definition: Objects.cpp:220
dd4hep::rec::MaterialManager::_pv
PlacedVolume _pv
Cached nodes.
Definition: MaterialManager.h:101
dd4hep::Material::A
double A() const
atomic number of the underlying material
Definition: Objects.cpp:187
dd4hep::rec::MaterialManager::_pos
Vector3D _pos
Definition: MaterialManager.h:104
dd4hep::Material::density
double density() const
density of the underlying material
Definition: Objects.cpp:198
dd4hep::rec::MaterialManager::~MaterialManager
~MaterialManager()
Definition: MaterialManager.cpp:31
dd4hep::Material
Handle class describing a material.
Definition: Objects.h:271
dd4hep::Volume
Handle class holding a placed volume (also called physical volume)
Definition: Volumes.h:371
dd4hep::rec::MaterialManager::_m
Material _m
Definition: MaterialManager.h:99
dd4hep::Material::Z
double Z() const
proton number of the underlying material
Definition: Objects.cpp:175
dd4hep::rec::MaterialManager::_p0
Vector3D _p0
cached last points
Definition: MaterialManager.h:104
dd4hep::rec::MaterialManager::_placeV
PlacementVec _placeV
Definition: MaterialManager.h:102
dd4hep::rec::MaterialManager::placementAt
PlacedVolume placementAt(const Vector3D &pos)
Definition: MaterialManager.cpp:193
dd4hep::rec::MaterialData
Definition: Material.h:33
dd4hep::rec::MaterialManager::entriesBetween
const ScanData entriesBetween(const Vector3D &p0, const Vector3D &p1, double eps=MaterialManager::epsilon)
As above, but optionally allow access to traversed placements.
Definition: MaterialManager.cpp:40
MaterialManager.h
dd4hep::rec::PlacementVec
std::vector< std::pair< PlacedVolume, double > > PlacementVec
Definition: MaterialManager.h:31
dd4hep::rec::MaterialManager::MaterialManager
MaterialManager()=delete
dd4hep
Namespace for the AIDA detector description toolkit.
Definition: AlignmentsCalib.h:28
dd4hep::rec::MaterialManager::_tgeoMgr
TGeoManager * _tgeoMgr
Reference to the TGeoManager.
Definition: MaterialManager.h:106
dd4hep::Material::radLength
double radLength() const
Access the radiation length of the underlying material.
Definition: Objects.cpp:209
MINSTEP
#define MINSTEP
Definition: MaterialManager.cpp:22
Exceptions.h