5 #include "TGeoVolume.h"
6 #include "TGeoManager.h"
8 #include "TVirtualGeoTrack.h"
27 if( ( p0 !=
_p0 ) || ( p1 !=
_p1 ) ) {
36 double startpoint[3], endpoint[3], direction[3];
38 for(
unsigned int i=0; i<3; i++) {
39 startpoint[i] = p0[i];
41 direction[i] = endpoint[i] - startpoint[i];
42 L+=direction[i]*direction[i];
44 double totDist = sqrt( L ) ;
47 for(
unsigned int i=0; i<3; i++)
48 direction[i]=direction[i]/totDist;
52 TGeoNode *node1 =
_tgeoMgr->InitTrack(startpoint, direction);
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.");
64 TGeoNode * node2 =
_tgeoMgr->FindNextBoundaryAndStep( 500, 1) ;
66 if( !node2 ||
_tgeoMgr->IsOutside() )
69 const double *position =
_tgeoMgr->GetCurrentPoint();
70 const double *previouspos =
_tgeoMgr->GetLastPoint();
74 TVirtualGeoTrack *track =
_tgeoMgr->GetLastTrack();
80 #if 1 //fg: is this still needed ?
84 position[1] +
MINSTEP * direction[1],
85 position[2] +
MINSTEP * direction[2] );
88 node2 =
_tgeoMgr->FindNextBoundaryAndStep(500, 1) ;
90 position =
_tgeoMgr->GetCurrentPoint();
91 previouspos =
_tgeoMgr->GetLastPoint();
99 double currDistance = ( posV - p0 ).r() ;
106 if( currDistance > totDist ) {
108 length = sqrt( pow(endpoint[0]-previouspos[0],2) +
109 pow(endpoint[1]-previouspos[1],2) +
110 pow(endpoint[2]-previouspos[2],2) );
112 track->AddPoint( endpoint[0], endpoint[1], endpoint[2], 0. );
115 if( length > epsilon )
116 _mV.push_back( std::make_pair(
Material( node1->GetMedium() ) , length ) ) ;
121 track->AddPoint( position[0], position[1], position[2], 0.);
123 if( length > epsilon )
124 _mV.push_back( std::make_pair(
Material( node1->GetMedium() ), length ) ) ;
132 _mV.push_back( std::make_pair(
Material( node1->GetMedium() ), totDist ) ) ;
154 TGeoNode *node=
_tgeoMgr->FindNode( pos[0], pos[1], pos[2] ) ;
157 std::stringstream err ;
158 err <<
" MaterialManager::material: No geometry node found at location: " << pos ;
159 throw std::runtime_error( err.str() );
174 std::stringstream sstr ;
177 double sum_rho_l = 0 ;
178 double sum_rho_l_over_A = 0 ;
179 double sum_rho_l_Z_over_A = 0 ;
181 double sum_l_over_x = 0 ;
183 double sum_l_over_lambda = 0 ;
185 for(
unsigned i=0,n=materials.size(); i<n ; ++i){
188 double l = materials[i].second ;
190 if( i != 0 ) sstr <<
"_" ;
191 sstr << mat.
name() <<
"_" << 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 ;
209 double rho = sum_rho_l / sum_l ;
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 ;
217 double x = sum_l / sum_l_over_x ;
220 double lambda = sum_l / sum_l_over_lambda ;
223 return MaterialData( sstr.str() , Z, A, rho, x, lambda ) ;
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.
const char * name() const
Access the object name (or "" if not supported by the object)
Handle class describing a material.
double intLength() const
Access the interaction length of the underlying material.
static LCDD & getInstance(void)
—Factory method----—
std::vector< std::pair< Material, double > > MaterialVec
DDSurfaces::Vector3D _pos
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.
double A() const
atomic number of the underlying material
double Z() const
proton number of the underlying material
double density() const
density of the underlying material