224{
225
226 static const GeoIntrusivePtr<const GeoBox> box1 = new GeoBox(1.*GeoModelKernelUnits::km,1*GeoModelKernelUnits::km,1*GeoModelKernelUnits::km);
227 static const GeoIntrusivePtr<const GeoShape>
s1 = [&] {
228 GeoTrf::Vector3D v1(0,0,-1*GeoModelKernelUnits::km);
229 GeoTrf::Transform3D ttt1 = GeoTrf::Transform3D::Identity()*GeoTrf::Translation3D(v1);
230 return new GeoShapeShift (&*box1, ttt1);
231 }();
232 static const GeoIntrusivePtr<const GeoShape>
s2 = [&] {
233 GeoTrf::Vector3D
v2(0,0,+1*GeoModelKernelUnits::km);
234 GeoTrf::Transform3D ttt2 = GeoTrf::Transform3D::Identity()*GeoTrf::Translation3D(v2);
235 return new GeoShapeShift (&*box1, ttt2);
236 }();
237
239 GeoTrf::Vector3D axis0(
v->GetPoint(0)-
v->GetPoint(1));
240 GeoTrf::Vector3D
axis(
v->GetPoint(1)-
v->GetPoint(0));
241 GeoTrf::Vector3D axis1 = GeoTrf::Vector3D::Identity();
242 GeoTrf::Vector3D axis2(
v->GetPoint(2)-
v->GetPoint(1));
244 double angle1=0;
245 double angle2=std::abs(std::atan2(
axis.cross(axis2).norm(),
axis.dot(axis2)))/2;
246 double delta_l1=0;
247 double delta_l2=
radius*std::tan(angle2);
248 double lengthnew=
length+delta_l2;
249 GeoShape*
solid=
new GeoTubs(0.,radius,lengthnew/2.,0.,2.0*
M_PI);
250
251 const GeoTrf::Vector3D
vt(0.,0.,-lengthnew/2.+delta_l2+2.);
252 GeoTrf::Transform3D rrr = GeoTrf::RotateY3D(angle2)*GeoTrf::RotateZ3D(
phi(axis2));
253 GeoShape *ssnew=new GeoShapeShift(&*s1,GeoTrf::Translation3D(vt)*rrr);
254
256
257 GeoTrf::Vector3D vref(0.,0.,-lengthnew/2.);
258 GeoTrf::Transform3D tref = GeoTrf::Transform3D::Identity()*GeoTrf::Translation3D(vref);
260 GeoTrf::Transform3D
r1 = GeoTrf::RotateZ3D(
phi(axis0))*GeoTrf::RotateY3D(
theta(axis0));
261 GeoTrf::Vector3D vtt(
v->GetPoint(0).x(),
v->GetPoint(0).y(),
v->GetPoint(0).z());
262 solid=
new GeoShapeShift(
solid,GeoTrf::Translation3D(vtt)*r1);
263
264 for (
int i=1;
i<
v->NrOfPoints()-1;
i++)
265 {
266 axis0=
v->GetPoint(i)-
v->GetPoint(i+1);
267 axis=
v->GetPoint(i+1)-
v->GetPoint(i);
268 axis1=
v->GetPoint(i)-
v->GetPoint(i-1);
269
271 angle1=std::abs(std::atan2(
axis.cross(axis1).norm(),
axis.dot(axis1)))/2;
272 delta_l1=
radius*std::tan(angle1);
273 delta_l2=0;
274 if (i<(
v->NrOfPoints()-2))
275 {
276 axis2=
v->GetPoint(i+2)-
v->GetPoint(i+1);
277 angle2=std::abs(std::atan2(
axis.cross(axis2).norm(),
axis.dot(axis2)))/2;
278 delta_l2=
radius*std::tan(angle2);
279 }
281 lengthnew=
length+delta_l1+delta_l2;
282
283 GeoTrf::Vector3D vvref(0.,0.,-lengthnew/2+delta_l1);
284 GeoTrf::Transform3D ttref = GeoTrf::Transform3D::Identity()*GeoTrf::Translation3D(vvref);
285
286 GeoShape*
ss=
new GeoTubs(0.,radius,lengthnew/2.,0.,2.0*
M_PI);
287
288 const GeoTrf::Vector3D vt1(0.,0.,+lengthnew/2.-delta_l1-2.);
289 const GeoTrf::Vector3D vt2(0.,0.,-lengthnew/2.+delta_l2+2.);
290 GeoTrf::Transform3D rrr1 = GeoTrf::RotateZ3D(-
phi(axis1))*GeoTrf::RotateY3D(angle1);
291 GeoTrf::Transform3D rrr2 = GeoTrf::RotateZ3D(
phi(axis2))*GeoTrf::RotateY3D(-angle2);
292 GeoTrf::Transform3D ttt1 = rrr1*GeoTrf::Translation3D(vt1);
293 GeoTrf::Transform3D ttt2 = rrr2*GeoTrf::Translation3D(vt2);
294 GeoShape *ssnew1=new GeoShapeShift(&*s2,ttt1);
295 ss =
new GeoShapeSubtraction(
ss,ssnew1);
296 if (i<(
v->NrOfPoints()-2))
297 {
298 GeoShape *ssnew2=new GeoShapeShift(&*s1,ttt2);
299 ss =
new GeoShapeSubtraction(
ss,ssnew2);
300 }
301
302 ss=
new GeoShapeShift(
ss,ttref);
303
304 GeoTrf::Transform3D rr1 = GeoTrf::RotateZ3D(
phi(axis0))*GeoTrf::RotateY3D(
theta(axis0));
305 const GeoTrf::Vector3D
vv(
v->GetPoint(i).x(),
v->GetPoint(i).y(),
v->GetPoint(i).z());
306 ss=
new GeoShapeShift(
ss,GeoTrf::Translation3D(vv)*rr1);
308 }
310}
Scalar phi() const
phi method
Scalar theta() const
theta method