DD4hep - The AIDA detector description toolkit for high energy physics experiments
DD4hep  Rev:Unversioneddirectory
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MaterialManager.cpp
Go to the documentation of this file.
2 #include "DD4hep/Exceptions.h"
3 #include "DD4hep/LCDD.h"
4 
5 #include "TGeoVolume.h"
6 #include "TGeoManager.h"
7 #include "TGeoNode.h"
8 #include "TVirtualGeoTrack.h"
9 
10 #define MINSTEP 1.e-5
11 
12 namespace DD4hep {
13  namespace DDRec {
14 
15 
16  MaterialManager::MaterialManager() : _mV(0), _m( Material() ), _p0(),_p1(),_pos() {
17 
18  _tgeoMgr = Geometry::LCDD::getInstance().world().volume()->GetGeoManager();
19  }
20 
22 
23  }
24 
26 
27  if( ( p0 != _p0 ) || ( p1 != _p1 ) ) {
28 
29  //---------------------------------------
30  _mV.clear() ;
31 
32  //
33  // algorithm copied from TGeoGearDistanceProperties.cc (A.Munnich):
34  //
35 
36  double startpoint[3], endpoint[3], direction[3];
37  double L=0;
38  for(unsigned int i=0; i<3; i++) {
39  startpoint[i] = p0[i];
40  endpoint[i] = p1[i];
41  direction[i] = endpoint[i] - startpoint[i];
42  L+=direction[i]*direction[i];
43  }
44  double totDist = sqrt( L ) ;
45 
46  //normalize direction
47  for(unsigned int i=0; i<3; i++)
48  direction[i]=direction[i]/totDist;
49 
50  _tgeoMgr->AddTrack(0, 12 ) ; // electron neutrino
51 
52  TGeoNode *node1 = _tgeoMgr->InitTrack(startpoint, direction);
53 
54  //check if there is a node at startpoint
55  if(!node1)
56  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.");
57 
58  while ( !_tgeoMgr->IsOutside() ) {
59 
60  // TGeoNode *node2;
61  // TVirtualGeoTrack *track;
62 
63  // step to (and over) the next Boundary
64  TGeoNode * node2 = _tgeoMgr->FindNextBoundaryAndStep( 500, 1) ;
65 
66  if( !node2 || _tgeoMgr->IsOutside() )
67  break;
68 
69  const double *position = _tgeoMgr->GetCurrentPoint();
70  const double *previouspos = _tgeoMgr->GetLastPoint();
71 
72  double length = _tgeoMgr->GetStep();
73 
74  TVirtualGeoTrack *track = _tgeoMgr->GetLastTrack();
75 
76  //protection against infinitive loop in root which should not happen, but well it does...
77  //work around until solution within root can be found when the step gets very small e.g. 1e-10
78  //and the next boundary is never reached
79 
80 #if 1 //fg: is this still needed ?
81  if( length < MINSTEP ) {
82 
83  _tgeoMgr->SetCurrentPoint( position[0] + MINSTEP * direction[0],
84  position[1] + MINSTEP * direction[1],
85  position[2] + MINSTEP * direction[2] );
86 
87  length = _tgeoMgr->GetStep();
88  node2 = _tgeoMgr->FindNextBoundaryAndStep(500, 1) ;
89 
90  position = _tgeoMgr->GetCurrentPoint();
91  previouspos = _tgeoMgr->GetLastPoint();
92  }
93 #endif
94  // printf( " -- step length : %1.8e %1.8e %1.8e %1.8e %1.8e %1.8e %1.8e - %s \n" , length ,
95  // position[0], position[1], position[2], previouspos[0], previouspos[1], previouspos[2] , node1->GetMedium()->GetMaterial()->GetName() ) ;
96 
97  DDSurfaces::Vector3D posV( position ) ;
98 
99  double currDistance = ( posV - p0 ).r() ;
100 
101  // //if the next boundary is further than end point
102  // if(fabs(position[0])>fabs(endpoint[0]) || fabs(position[1])>fabs(endpoint[1])
103  // || fabs(position[2])>fabs(endpoint[2]))
104 
105  //if we travelled too far:
106  if( currDistance > totDist ) {
107 
108  length = sqrt( pow(endpoint[0]-previouspos[0],2) +
109  pow(endpoint[1]-previouspos[1],2) +
110  pow(endpoint[2]-previouspos[2],2) );
111 
112  track->AddPoint( endpoint[0], endpoint[1], endpoint[2], 0. );
113 
114 
115  if( length > epsilon )
116  _mV.push_back( std::make_pair( Material( node1->GetMedium() ) , length ) ) ;
117 
118  break;
119  }
120 
121  track->AddPoint( position[0], position[1], position[2], 0.);
122 
123  if( length > epsilon )
124  _mV.push_back( std::make_pair( Material( node1->GetMedium() ), length ) ) ;
125 
126  node1 = node2 ;
127  }
128 
129 
130  //fg: protect against empty list:
131  if( _mV.empty() ){
132  _mV.push_back( std::make_pair( Material( node1->GetMedium() ), totDist ) ) ;
133  }
134 
135 
136  _tgeoMgr->ClearTracks();
137 
138  _tgeoMgr->CleanGarbage();
139 
140  //---------------------------------------
141 
142  _p0 = p0 ;
143  _p1 = p1 ;
144  }
145 
146  return _mV ; ;
147  }
148 
149 
151 
152  if( pos != _pos ) {
153 
154  TGeoNode *node=_tgeoMgr->FindNode( pos[0], pos[1], pos[2] ) ;
155 
156  if( ! node ) {
157  std::stringstream err ;
158  err << " MaterialManager::material: No geometry node found at location: " << pos ;
159  throw std::runtime_error( err.str() );
160  }
161 
162  // std::cout << " @@@ MaterialManager::material @ " << pos << " found volume : " << node->GetName() << std::endl ;
163 
164  _m = Material( node->GetMedium() ) ;
165 
166  _pos = pos ;
167  }
168 
169  return _m ; ;
170  }
171 
173 
174  std::stringstream sstr ;
175 
176  double sum_l = 0 ;
177  double sum_rho_l = 0 ;
178  double sum_rho_l_over_A = 0 ;
179  double sum_rho_l_Z_over_A = 0 ;
180  //double sum_rho_l_over_x = 0 ;
181  double sum_l_over_x = 0 ;
182  //double sum_rho_l_over_lambda = 0 ;
183  double sum_l_over_lambda = 0 ;
184 
185  for(unsigned i=0,n=materials.size(); i<n ; ++i){
186 
187  Material mat = materials[i].first ;
188  double l = materials[i].second ;
189 
190  if( i != 0 ) sstr << "_" ;
191  sstr << mat.name() << "_" << l ;
192 
193  double rho = mat.density() ;
194  double A = mat.A() ;
195  double Z = mat.Z() ;
196  double x = mat.radLength() ;
197  double lambda = mat.intLength() ;
198 
199  sum_l += l ;
200  sum_rho_l += rho * l ;
201  sum_rho_l_over_A += rho * l / A ;
202  sum_rho_l_Z_over_A += rho * l * Z / A ;
203  sum_l_over_x += l / x ;
204  sum_l_over_lambda += l / lambda ;
205  // sum_rho_l_over_x += rho * l / x ;
206  // sum_rho_l_over_lambda += rho * l / lambda ;
207  }
208 
209  double rho = sum_rho_l / sum_l ;
210 
211  double A = sum_rho_l / sum_rho_l_over_A ;
212  double Z = sum_rho_l_Z_over_A / sum_rho_l_over_A ;
213 
214  // radiation and interaction lengths already given in cm - average by length
215 
216  // double x = sum_rho_l / sum_rho_l_over_x ;
217  double x = sum_l / sum_l_over_x ;
218 
219  // double lambda = sum_rho_l / sum_rho_l_over_lambda ;
220  double lambda = sum_l / sum_l_over_lambda ;
221 
222 
223  return MaterialData( sstr.str() , Z, A, rho, x, lambda ) ;
224 
225  }
226 
227  } /* namespace DDRec */
228 } /* namespace DD4hep */
virtual DetElement world() const =0
Return reference to the top-most (world) detector element.
Volume volume() const
Access to the logical volume of the detector element's placement.
Definition: Detector.cpp:311
const char * name() const
Access the object name (or "" if not supported by the object)
Definition: Handle.inl:36
Handle class describing a material.
Definition: Objects.h:300
double intLength() const
Access the interaction length of the underlying material.
Definition: Objects.cpp:225
#define MINSTEP
static LCDD & getInstance(void)
—Factory method----—
Definition: LCDDImp.cpp:87
std::vector< std::pair< Material, double > > MaterialVec
const Material & materialAt(const DDSurfaces::Vector3D &pos)
const MaterialVec & materialsBetween(const DDSurfaces::Vector3D &p0, const DDSurfaces::Vector3D &p1, double epsilon=1e-4)
MaterialData createAveragedMaterial(const MaterialVec &materials)
double radLength() const
Access the radiation length of the underlying material.
Definition: Objects.cpp:213
double A() const
atomic number of the underlying material
Definition: Objects.cpp:189
double Z() const
proton number of the underlying material
Definition: Objects.cpp:178
double density() const
density of the underlying material
Definition: Objects.cpp:201