153{
154 G4VSolid* theSolid(nullptr);
155
158 if(sharedShapes.find(geoShape)!=sharedShapes.end())
159 return sharedShapes[geoShape];
160
161
162 G4VSolid* solidA(nullptr);
163 G4VSolid* solidB(nullptr);
165
166
167
168
169 if(geoShape->typeID() == GeoBox::getClassTypeID() )
170 {
171 const GeoBox* theBox = dynamic_cast<const GeoBox*> (geoShape);
172 if (nullptr==theBox) throw std::runtime_error("TypeID did not match cast for box");
173 if (
n.empty())
n=
"G4Box";
174 if (theBox->getXHalfLength()<=0.){
ATH_MSG_WARNING(
"Box " << n <<
" has an x side of " << theBox->getXHalfLength() <<
" - using std::abs.");}
175 if (theBox->getYHalfLength()<=0.){
ATH_MSG_WARNING(
"Box " << n <<
" has an y side of " << theBox->getYHalfLength() <<
" - using std::abs.");}
176 if (theBox->getZHalfLength()<=0.){
ATH_MSG_WARNING(
"Box " << n <<
" has an z side of " << theBox->getZHalfLength() <<
" - using std::abs.");}
177 theSolid = new G4Box(n,
178 std::abs(theBox->getXHalfLength()),
179 std::abs(theBox->getYHalfLength()),
180 std::abs(theBox->getZHalfLength()));
181 }
182
183
184
185 else if(geoShape->typeID() == GeoTube::getClassTypeID() )
186 {
187 const GeoTube* theTube = dynamic_cast<const GeoTube*> (geoShape);
188 if (nullptr==theTube) throw std::runtime_error("TypeID did not match cast for tube");
189 if (
n.empty())
n=
"G4Tube";
190 if (theTube->getRMax()<=0.){
ATH_MSG_WARNING(
"Tube " << n <<
" has a max radius of " << theTube->getRMax() <<
" - using std::abs.");}
191 if (theTube->getZHalfLength()<=0.){
ATH_MSG_WARNING(
"Tube " << n <<
" has a z half length of " << theTube->getZHalfLength() <<
" - using std::abs.");}
192 if (theTube->getRMax()<theTube->getRMin()){
ATH_MSG_WARNING(
"Tube " << n <<
" has a max radius of " << theTube->getRMax() <<
" and a min radius of " << theTube->getRMin());}
193 theSolid = new G4Tubs(n,
194 theTube->getRMin(),
195 std::abs(theTube->getRMax()),
196 std::abs(theTube->getZHalfLength()),
197 0.,360*CLHEP::deg);
198 }
199
200
201
202 else if(geoShape->typeID() == GeoTubs::getClassTypeID() )
203 {
204 const GeoTubs* theTubs = dynamic_cast<const GeoTubs*> (geoShape);
205 if (nullptr==theTubs) throw std::runtime_error("TypeID did not match cast for tubs");
206 if (
n.empty())
n=
"G4Tubs";
207 if (theTubs->getRMin()<0.){
ATH_MSG_WARNING(
"Tubs " << n <<
" has a min radius of " << theTubs->getRMax());}
208 if (theTubs->getRMax()<=0.){
ATH_MSG_WARNING(
"Tubs " << n <<
" has a max radius of " << theTubs->getRMax() <<
" - using std::abs.");}
209 if (theTubs->getZHalfLength()<=0.){
ATH_MSG_WARNING(
"Tubs " << n <<
" has a half length of " << theTubs->getZHalfLength() <<
" - using std::abs.");}
210 if (theTubs->getRMax()<theTubs->getRMin()){
ATH_MSG_WARNING(
"Tubs " << n <<
" has a max radius of " << theTubs->getRMax() <<
" and a min radius of " << theTubs->getRMin());}
211 if (theTubs->getDPhi()<=0.){
ATH_MSG_WARNING(
"Tubs " << n <<
" has a dPhi of " << theTubs->getDPhi());}
212 theSolid = new G4Tubs(n,
213 theTubs->getRMin(),
214 std::abs(theTubs->getRMax()),
215 std::abs(theTubs->getZHalfLength()),
216 theTubs->getSPhi(),
217 theTubs->getDPhi());
218 }
219
220
221
222 else if(geoShape->typeID() == GeoTrd::getClassTypeID() )
223 {
224 const GeoTrd* theTrd = dynamic_cast<const GeoTrd*> (geoShape);
225 if (nullptr==theTrd) throw std::runtime_error("TypeID did not match cast for trd");
226 if (
n.empty())
n=
"G4Trd";
227 if (theTrd->getXHalfLength1()<0.){
ATH_MSG_WARNING(
"Trd " << n <<
" has a x half length 1 of " << theTrd->getXHalfLength1() <<
" - using std::abs.");}
228 if (theTrd->getXHalfLength2()<0.){
ATH_MSG_WARNING(
"Trd " << n <<
" has a x half length 2 of " << theTrd->getXHalfLength2() <<
" - using std::abs.");}
229 if (theTrd->getYHalfLength1()<0.){
ATH_MSG_WARNING(
"Trd " << n <<
" has a y half length 1 of " << theTrd->getYHalfLength1() <<
" - using std::abs.");}
230 if (theTrd->getYHalfLength2()<0.){
ATH_MSG_WARNING(
"Trd " << n <<
" has a y half length 2 of " << theTrd->getYHalfLength2() <<
" - using std::abs.");}
231 if (theTrd->getZHalfLength()<=0.){
ATH_MSG_WARNING(
"Trd " << n <<
" has a z half length of " << theTrd->getZHalfLength() <<
" - using std::abs.");}
232 if (theTrd->getXHalfLength1()<=0. && theTrd->getXHalfLength2()<=0.){
ATH_MSG_WARNING(
"Trd " << n <<
" has an x half length 1 of " << theTrd->getXHalfLength1()
233 << " and an x half length 2 of " << theTrd->getXHalfLength2() << " - using std::abs.");}
234 if (theTrd->getYHalfLength1()<=0. && theTrd->getYHalfLength2()<=0.){
ATH_MSG_WARNING(
"Trd " << n <<
" has a y half length 1 of " << theTrd->getYHalfLength1()
235 << " and a y half length 2 of " << theTrd->getYHalfLength2() << " - using std::abs.");}
236 theSolid = new G4Trd(n,
237 std::abs(theTrd->getXHalfLength1()),
238 std::abs(theTrd->getXHalfLength2()),
239 std::abs(theTrd->getYHalfLength1()),
240 std::abs(theTrd->getYHalfLength2()),
241 std::abs(theTrd->getZHalfLength()));
242 }
243
244
245
246 else if(geoShape->typeID() == GeoPcon::getClassTypeID())
247 {
248 const GeoPcon* thePcon = dynamic_cast<const GeoPcon*>(geoShape);
249 if (nullptr==thePcon) throw std::runtime_error("TypeID did not match cast for pcon");
250 theSolid =
new G4Polycone(
n.empty()?
"G4Polycone":std::move(n),
251 thePcon->getSPhi(),
252 thePcon->getDPhi(),
253 thePcon->getNPlanes(),
254 thePcon->getZPlaneBuff(),
255 thePcon->getRMinBuff(),
256 thePcon->getRMaxBuff());
257 }
258
259
260
261 else if(geoShape->typeID() == GeoCons::getClassTypeID())
262 {
263 const GeoCons* theCons = dynamic_cast<const GeoCons*>(geoShape);
264 if (nullptr==theCons) throw std::runtime_error("TypeID did not match cast for cons");
265 if (
n.empty())
n=
"G4Cons";
266 if (theCons->getRMax1()<0.){
ATH_MSG_WARNING(
"Cons " << n <<
" has a max radius 1 of " << theCons->getRMax1() <<
" - will use std::abs.");}
267 if (theCons->getRMax2()<0.){
ATH_MSG_WARNING(
"Cons " << n <<
" has a max radius 2 of " << theCons->getRMax2() <<
" - will use std::abs.");}
268 if (theCons->getRMin1()<0.){
ATH_MSG_WARNING(
"Cons " << n <<
" has a min radius 1 of " << theCons->getRMin1());}
269 if (theCons->getRMin2()<0.){
ATH_MSG_WARNING(
"Cons " << n <<
" has a min radius 2 of " << theCons->getRMin2());}
270 if (theCons->getDZ()<=0){
ATH_MSG_WARNING(
"Cons " << n <<
" has a DZ of " << theCons->getDZ() <<
" - will use std::abs.");}
271 if (theCons->getRMax1()<=0. && theCons->getRMax2()<=0.){
ATH_MSG_WARNING(
"Cons " << n <<
" has a max radius 1 of " << theCons->getRMax1()
272 << " and a max radius 2 of " << theCons->getRMax2() << " - will use std::abs.");}
273 theSolid = new G4Cons(n,
274 theCons->getRMin1(),
275 std::abs(theCons->getRMax1()),
276 theCons->getRMin2(),
277 std::abs(theCons->getRMax2()),
278 std::abs(theCons->getDZ()),
279 theCons->getSPhi(),
280 theCons->getDPhi());
281 }
282
283
284
285 else if(geoShape->typeID() == GeoPara::getClassTypeID())
286 {
287 const GeoPara* thePara = dynamic_cast<const GeoPara*>(geoShape);
288 if (nullptr==thePara) throw std::runtime_error("TypeID did not match cast for para");
289 if (
n.empty())
n=
"G4Para";
290 if (thePara->getXHalfLength()<=0.){
ATH_MSG_WARNING(
"Para " << n <<
" has an x side of " << thePara->getXHalfLength() <<
" - using std::abs.");}
291 if (thePara->getYHalfLength()<=0.){
ATH_MSG_WARNING(
"Para " << n <<
" has an y side of " << thePara->getYHalfLength() <<
" - using std::abs.");}
292 if (thePara->getZHalfLength()<=0.){
ATH_MSG_WARNING(
"Para " << n <<
" has an z side of " << thePara->getZHalfLength() <<
" - using std::abs.");}
293 theSolid = new G4Para(n,
294 std::abs(thePara->getXHalfLength()),
295 std::abs(thePara->getYHalfLength()),
296 std::abs(thePara->getZHalfLength()),
297 thePara->getAlpha(),
298 thePara->getTheta(),
299 thePara->getPhi());
300 }
301
302
303
304 else if(geoShape->typeID() == GeoPgon::getClassTypeID())
305 {
306 const GeoPgon* thePgon = dynamic_cast<const GeoPgon*>(geoShape);
307 if (nullptr==thePgon) throw std::runtime_error("TypeID did not match cast for pgon");
308
309
310 unsigned nPlanes = thePgon->getNPlanes();
311 auto rInner = std::make_unique<double[]>(nPlanes);
312 auto rOuter = std::make_unique<double[]>(nPlanes);
313 double alpha = thePgon->getDPhi()/(2*thePgon->getNSides());
317 }
318
319 theSolid =
new G4Polyhedra(
n.empty()?
"G4Polyhedra":std::move(n),
320 thePgon->getSPhi(),
321 thePgon->getDPhi(),
322 thePgon->getNSides(),
323 nPlanes,
324 thePgon->getZPlaneBuff(),
325 rInner.get(),
326 rOuter.get());
327 }
328
329
330
331 else if(geoShape->typeID() == GeoTrap::getClassTypeID())
332 {
333 const GeoTrap* theTrap = dynamic_cast<const GeoTrap*>(geoShape);
334 if (nullptr==theTrap) throw std::runtime_error("TypeID did not match cast for trap");
335 if (
n.empty())
n=
"G4Trap";
336 if (theTrap->getZHalfLength()<=0.){
ATH_MSG_WARNING(
"Trap " << n <<
" has an z side of " << theTrap->getZHalfLength() <<
" - using std::abs.");}
337 theSolid = new G4Trap(n,
338 std::abs(theTrap->getZHalfLength()),
339 theTrap->getTheta(),
340 theTrap->getPhi(),
341 theTrap->getDydzn(),
342 theTrap->getDxdyndzn(),
343 theTrap->getDxdypdzn(),
344 theTrap->getAngleydzn(),
345 theTrap->getDydzp(),
346 theTrap->getDxdyndzp(),
347 theTrap->getDxdypdzp(),
348 theTrap->getAngleydzp());
349 }
350
351
352
353 else if(geoShape->typeID() == GeoSimplePolygonBrep::getClassTypeID())
354 {
355 const GeoSimplePolygonBrep* theBrep = dynamic_cast<const GeoSimplePolygonBrep*>(geoShape);
356 if (nullptr==theBrep) throw std::runtime_error("TypeID did not match cast for brep");
357 if (
n.empty())
n=
"G4ExtrudedSolid";
358 double dz = theBrep->getDZ();
359 int nVertices = theBrep->getNVertices();
360
361 G4TwoVector off(0,0);
362 std::vector<G4TwoVector> polygon;
363
364 polygon.reserve(nVertices);
365for(
int i=0;
i<nVertices;
i++)
366 polygon.push_back(G4TwoVector(theBrep->getXVertex(nVertices-1-i),theBrep->getYVertex(nVertices-1-i)));
367
368 theSolid = new G4ExtrudedSolid(n,polygon,dz,off,1,off,1);
369 }
370
371
372
373 else if(geoShape->typeID() == GeoTessellatedSolid::getClassTypeID())
374 {
375 const GeoTessellatedSolid* theTessellated = dynamic_cast<const GeoTessellatedSolid*>(geoShape);
376 if (nullptr==theTessellated) throw std::runtime_error("TypeID did not match cast for tessellated solid");
377 if(
n.empty())
n=
"G4TessellatedSolid";
378
379 G4TessellatedSolid* g4Tessellated = new G4TessellatedSolid(n);
380 for(
size_t i=0;
i<theTessellated->getNumberOfFacets(); ++
i) {
381 GeoFacet* geoFacet = theTessellated->getFacet(i);
382 G4FacetVertexType
vertexType = (geoFacet->getVertexType()==GeoFacet::ABSOLUTE? ABSOLUTE : RELATIVE);
383 G4VFacet* g4Facet(nullptr);
384 if(geoFacet->getNumberOfVertices()==3)
388 vertexType);
389 else
394 vertexType);
395
396 g4Tessellated->AddFacet(g4Facet);
397 }
398 g4Tessellated->SetSolidClosed(true);
399 theSolid = g4Tessellated;
400 }
401
402
403
404 else if(geoShape->typeID() == GeoEllipticalTube::getClassTypeID())
405 {
406 const GeoEllipticalTube* theEltube = dynamic_cast<const GeoEllipticalTube*>(geoShape);
407 if (nullptr==theEltube) throw std::runtime_error("TypeID did not match cast for elliptical tube");
408 if (
n.empty())
n=
"G4EllipticalTube";
409
410 if (theEltube->getXHalfLength()<=0.){
ATH_MSG_WARNING(
"Eltube " << n <<
" has an x side of " << theEltube->getXHalfLength() <<
" - using std::abs.");}
411 if (theEltube->getYHalfLength()<=0.){
ATH_MSG_WARNING(
"Eltube " << n <<
" has an y side of " << theEltube->getYHalfLength() <<
" - using std::abs.");}
412 if (theEltube->getZHalfLength()<=0.){
ATH_MSG_WARNING(
"Eltube " << n <<
" has an z side of " << theEltube->getZHalfLength() <<
" - using std::abs.");}
413 G4EllipticalTube* g4Eltube = new G4EllipticalTube(n
414 ,std::abs(theEltube->getXHalfLength())
415 ,std::abs(theEltube->getYHalfLength())
416 ,std::abs(theEltube->getZHalfLength()));
417 theSolid = g4Eltube;
418 }
419
420
421
422 else if(geoShape->typeID() == GeoTorus::getClassTypeID() ) {
423 const GeoTorus* theTorus = dynamic_cast<const GeoTorus*> (geoShape);
424 if (nullptr==theTorus) throw std::runtime_error("TypeID did not match cast for torus");
425 if (
n.empty())
n=
"G4Torus";
426
427 theSolid = new G4Torus(n,
428 theTorus->getRMin(),
429 theTorus->getRMax(),
430 theTorus->getRTor(),
431 theTorus->getSPhi(),
432 theTorus->getDPhi());
433 }
434
435
436
437 else if(geoShape->typeID() == GeoGenericTrap::getClassTypeID()) {
438 const GeoGenericTrap* theGenTrap = dynamic_cast<const GeoGenericTrap*>(geoShape);
439 if (nullptr==theGenTrap) throw std::runtime_error("TypeID did not match cast for generic trap");
440 if (
n.empty())
n=
"G4GenericTrap";
441 if (theGenTrap->getZHalfLength()<=0.){
ATH_MSG_WARNING(
"GenTrap " << n <<
" has an z side of " << theGenTrap->getZHalfLength() <<
" - using std::abs.");}
442
443
444 std::vector<CLHEP::Hep2Vector> clhepVertices;
445 clhepVertices.reserve(theGenTrap->getVertices().size());
446 for(const GeoTrf::Vector2D& geoVertex : theGenTrap->getVertices()) {
447 clhepVertices.push_back(CLHEP::Hep2Vector(geoVertex.x(),geoVertex.y()));
448 }
449
450 G4GenericTrap* g4GenTrap = new G4GenericTrap(n
451 ,std::abs(theGenTrap->getZHalfLength())
452 ,clhepVertices);
453 theSolid = g4GenTrap;
454 }
455
456
457
458
459
460
461
462 else if (geoShape->typeID() == GeoShapeShift::getClassTypeID() )
463 {
464 const GeoShapeShift* theShapeShift = dynamic_cast<const GeoShapeShift*> (geoShape);
465 if (nullptr==theShapeShift) throw std::runtime_error("TypeID did not match cast for shape shift");
466 if (
n.empty())
n=
"DisplacedSolid";
467 G4VSolid * undisplacedSolid = Build(theShapeShift->getOp());
469 }
470
471
472
473 else if (geoShape->typeID() == GeoShapeUnion::getClassTypeID() )
474 {
475 const GeoShapeUnion* theUnion = dynamic_cast<const GeoShapeUnion*> (geoShape);
476 if (nullptr==theUnion) throw std::runtime_error("TypeID did not match cast for union");
477 if (
n.empty())
n=
"Union";
478 solidA = Build(theUnion->getOpA());
479 solidB = Build(theUnion->getOpB());
480 theSolid = new G4UnionSolid(n, solidA, solidB);
481 }
482
483
484
485 else if (geoShape->typeID() == GeoShapeIntersection::getClassTypeID() )
486 {
487 const GeoShapeIntersection* theIntersection = dynamic_cast<const GeoShapeIntersection*>(geoShape);
488 if (nullptr==theIntersection) throw std::runtime_error("TypeID did not match cast for intersection");
489 if (
n.empty())
n=
"Intersection";
490 solidA = Build(theIntersection->getOpA());
491 solidB = Build(theIntersection->getOpB());
492 theSolid = new G4IntersectionSolid(n, solidA, solidB);
493 }
494
495
496
497 else if (geoShape->typeID() == GeoShapeSubtraction::getClassTypeID() )
498 {
499 const GeoShapeSubtraction* theSubtraction = dynamic_cast<const GeoShapeSubtraction*>(geoShape);
500 if (nullptr==theSubtraction) throw std::runtime_error("TypeID did not match cast for subtraction");
501 if (
n.empty())
n=
"Subtraction";
502 solidA = Build(theSubtraction->getOpA());
503 solidB = Build(theSubtraction->getOpB());
504 theSolid = new G4SubtractionSolid(n, solidA, solidB);
505 }
506
507
508
509 else if(geoShape->typeID() == GeoUnidentifiedShape::getClassTypeID())
510 {
511 const GeoUnidentifiedShape* customShape = dynamic_cast<const GeoUnidentifiedShape*> (geoShape);
512 if (nullptr==customShape) throw std::runtime_error("TypeID did not match cast for custom shape");
513 if (customShape->name()=="LArCustomShape") {
514
515
516
517 ISvcLocator* svcLocator = Gaudi::svcLocator();
518
519 SmartIF<IGeoModelSvc> geoModel{svcLocator->service ("GeoModelSvc")};
520 if ( !geoModel ) {
521 G4Exception(
522 "Geo2G4SolidFactory", "AccessGeoModel", FatalException,
523 "Build cannot access GeoModelSvc");
524 }
525 SmartIF<IGeoDbTagSvc> geoDbTagSvc{svcLocator->service ("GeoDbTagSvc")};
526 if ( !geoDbTagSvc ) {
527 G4Exception(
528 "Geo2G4SolidFactory", "AccessDbTagSvc", FatalException,
529 "Build cannot access DbTagSvc");
530 }
531
532 SmartIF<IRDBAccessSvc> pAccessSvc{svcLocator->service(geoDbTagSvc->getParamSvcName())};
533 if ( !pAccessSvc ) {
534 G4Exception(
535 "Geo2G4SolidFactory", "AccessAccessSvc", FatalException,
536 "Build cannot access AccessSvc");
537 }
540
541
542
543 std::string customName = customShape->asciiData();
544 customSolidMap::const_iterator
it = customSolids.find(customName);
545 if(it!=customSolids.end())
546 theSolid =
it->second;
547 else
548 {
549 theSolid = nullptr;
550 if(customName.find("Slice") != std::string::npos){
551 theSolid = createLArWheelSliceSolid(customShape,emecData);
552 } else {
553 theSolid = createLArWheelSolid(customName, s_lwsTypes.at(customName) ,emecData);
554 }
555 if ( nullptr == theSolid ) {
556 std::string
error = std::string(
"Can't create LArWheelSolid for name ") + customName +
" in Geo2G4SolidFactory::Build";
557 throw std::runtime_error(
error);
558 }
559
560 if(theSolid != nullptr) customSolids[customName] = theSolid;
561 }
562 }
563 }
564
565
566
567 else
568 {
569 ATH_MSG_FATAL(
"Sorry this solid is not yet implemented... ");
572 return nullptr;
573 }
574
575 sharedShapes[geoShape] = theSolid;
576 return theSolid;
577}
#define ATH_MSG_WARNING(x)
std::map< const GeoShape *, G4VSolid *, std::less< const GeoShape * > > shapesMap
std::map< std::string, G4VSolid *, std::less< std::string > > customSolidMap
This is a helper class to query the version tags from GeoModelSvc and determine the appropriate tag a...
HepGeom::Transform3D EigenTransformToCLHEP(const Amg::Transform3D &eigenTransf)
Converts an Eigen-based Amg::Transform3D into a CLHEP-based HepGeom::Transform3D.
CLHEP::Hep3Vector EigenToHep3Vector(const Amg::Vector3D &eigenvector)
Converts an Eigen-based Amg::Vector3D into a CLHEP-based CLHEP::Hep3Vector.
EMECData toEMECData(IRDBAccessSvc *rdbAccess, const DecodeVersionKey &larVersionKey)