Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
#include "DDRec/MaterialManager.h"
#include "DD4hep/Exceptions.h"
#include "DD4hep/LCDD.h"
#include "TGeoVolume.h"
#include "TGeoManager.h"
#include "TGeoNode.h"
#include "TVirtualGeoTrack.h"
#define MINSTEP 1.e-5
namespace DD4hep {
namespace DDRec {
MaterialManager::MaterialManager() : _mV(0), _m( Material() ), _p0(),_p1(),_pos() {
_tgeoMgr = Geometry::LCDD::getInstance().world().volume()->GetGeoManager();
}
MaterialManager::~MaterialManager(){
}
const MaterialVec&MaterialManager:: materials(const DDSurfaces::Vector3D& p0, const DDSurfaces::Vector3D& p1 ) {
if( ( p0 != _p0 ) && ( p1 != _p1 ) ) {
//---------------------------------------
_mV.clear() ;
//
// algorithm copied from TGeoGearDistanceProperties.cc :
//
double startpoint[3], endpoint[3], direction[3];
double L=0;
for(unsigned int i=0; i<3; i++) {
startpoint[i] = p0[i];
endpoint[i] = p1[i];
direction[i] = endpoint[i] - startpoint[i];
L+=direction[i]*direction[i];
}
double totDist = sqrt( L ) ;
//normalize direction
for(unsigned int i=0; i<3; i++)
direction[i]=direction[i]/sqrt(L);
_tgeoMgr->AddTrack(0,12);
TGeoNode *node1 = _tgeoMgr->InitTrack(startpoint, direction);
//check if there is a node at startpoint
if(!node1)
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.");
while ( !_tgeoMgr->IsOutside() ) {
TGeoNode *node2;
TVirtualGeoTrack *track;
node2 = _tgeoMgr->FindNextBoundaryAndStep(500, 1) ;
if(!node2 || _tgeoMgr->IsOutside())
break;
const double *position =(double*) _tgeoMgr->GetCurrentPoint();
const double *previouspos =(double*) _tgeoMgr->GetLastPoint();
double length=_tgeoMgr->GetStep();
track = _tgeoMgr->GetLastTrack();
//protection against infinitive loop in root which should not happen, but well it does...
//work around until solution within root can be found when the step gets very small e.g. 1e-10
//and the next boundary is never reached
if(length<MINSTEP) {
_tgeoMgr->SetCurrentPoint(position[0]+MINSTEP*direction[0], position[1]+MINSTEP*direction[1], position[2]+MINSTEP*direction[2]);
length=_tgeoMgr->GetStep();
node2 = _tgeoMgr->FindNextBoundaryAndStep(500, 1) ;
//fg: need to update positions as well !?
position = (const double*) _tgeoMgr->GetCurrentPoint();
previouspos = (const double*) _tgeoMgr->GetLastPoint();
}
// printf( " -- step length : %1.8e %1.8e %1.8e %1.8e %1.8e %1.8e %1.8e - %s \n" , length ,
// position[0], position[1], position[2], previouspos[0], previouspos[1], previouspos[2] , node1->GetMedium()->GetMaterial()->GetName() ) ;
DDSurfaces::Vector3D posV( position ) ;
double currDistance = ( posV - p0 ).r() ;
// //if the next boundary is further than end point
// if(fabs(position[0])>fabs(endpoint[0]) || fabs(position[1])>fabs(endpoint[1])
// || fabs(position[2])>fabs(endpoint[2]))
//if we travelled too far:
if( currDistance > totDist ) {
length=sqrt(pow(endpoint[0]-previouspos[0],2)
+pow(endpoint[1]-previouspos[1],2)
+pow(endpoint[2]-previouspos[2],2));
track->AddPoint(endpoint[0], endpoint[1], endpoint[2], 0.);
_mV.push_back( std::make_pair( Material( node1->GetMedium() ) , length ) ) ;
break;
}
track->AddPoint(position[0], position[1], position[2], 0.);
_mV.push_back( std::make_pair( Material( node1->GetMedium() ), length ) ) ;
node1=node2;
}
_tgeoMgr->ClearTracks();
_tgeoMgr->CleanGarbage();
//---------------------------------------
_p0 = p0 ;
_p1 = p1 ;
}
return _mV ; ;
}
const Geometry::Material& MaterialManager::material(const DDSurfaces::Vector3D& pos ){
if( pos != _pos ) {
TGeoNode *node=_tgeoMgr->FindNode( pos[0], pos[1], pos[2] ) ;
if( ! node ) {
std::stringstream err ;
err << " MaterialManager::material: No geometry node found at location: " << pos ;
throw std::runtime_error( err.str() );
}
// std::cout << " @@@ MaterialManager::material @ " << pos << " found volume : " << node->GetName() << std::endl ;
_m = Material( node->GetMedium() ) ;
_pos = pos ;
}
return _m ; ;
}
} /* namespace DDRec */
} /* namespace DD4hep */