u4/tests/TubePhiCutQuarterShellTest.cc

0001 #include <cassert>
0002 #include <cmath>
0003 #include <cstdlib>
0004 #include <iostream>
0005 #include <vector>
0006
0007 #include <plog/Appenders/ConsoleAppender.h>
0008 #include <plog/Formatters/TxtFormatter.h>
0009 #include <plog/Init.h>
0010 #include <plog/Log.h>
0011
0012 using plog::error;
0013 using plog::fatal;
0014 using plog::info;
0015
0016 #include "G4LogicalVolume.hh"
0017 #include "G4PhysicalConstants.hh"
0018 #include "G4SystemOfUnits.hh"
0019 #include "G4Tubs.hh"
0020 #include "G4VPhysicalVolume.hh"
0021
0022 #include "config_path.h"
0023
0024 #include "U4GDML.h"
0025 #include "U4Solid.h"
0026 #include "U4Volume.h"
0027
0028 #include "s_csg.h"
0029
0030 namespace
0031 {
0032 inline constexpr char testGeomFile[] = SIMPHONY_TEST_GEOM_DIR "/tube_phicut_quarter_shell.gdml";
0033
0034 bool close(double actual, double expected, double tolerance = 1.e-9)
0035 {
0036     return std::fabs(actual - expected) <= tolerance * (1. + std::fabs(expected));
0037 }
0038 } // namespace
0039
0040 int main()
0041 {
0042     static plog::ConsoleAppender<plog::TxtFormatter> consoleAppender;
0043     plog::init(plog::info, &consoleAppender);
0044
0045     const G4VPhysicalVolume* world = U4GDML::Read(testGeomFile);
0046     if (world == nullptr)
0047     {
0048         std::cerr << "failed to load " << testGeomFile << std::endl;
0049         return EXIT_FAILURE;
0050     }
0051
0052     const G4VPhysicalVolume* tubePV = U4Volume::FindPV(world, "QuarterTube_pv");
0053     const G4LogicalVolume*   tubeLV = tubePV ? tubePV->GetLogicalVolume() : nullptr;
0054     const G4Tubs*            tube = tubeLV ? dynamic_cast<const G4Tubs*>(tubeLV->GetSolid()) : nullptr;
0055     if (tube == nullptr)
0056     {
0057         std::cerr << "failed to find partial-phi G4Tubs QuarterTube_pv" << std::endl;
0058         return EXIT_FAILURE;
0059     }
0060
0061     const double startPhi = tube->GetStartPhiAngle() / CLHEP::radian;
0062     const double deltaPhi = tube->GetDeltaPhiAngle() / CLHEP::radian;
0063     if (!close(startPhi, 0.) || !close(deltaPhi, 0.5 * CLHEP::pi))
0064     {
0065         std::cerr << "GDML tube did not preserve its quarter-phi interval" << std::endl;
0066         return EXIT_FAILURE;
0067     }
0068
0069     s_csg* csg = new s_csg;
0070     assert(csg);
0071
0072     const int lvid = 0;
0073     sn*       root = U4Solid::Convert(tube, lvid, 0, 1);
0074     if (root == nullptr || root->typecode != CSG_INTERSECTION)
0075     {
0076         std::cerr << "partial annular G4Tubs did not convert to a canonical intersection" << std::endl;
0077         return EXIT_FAILURE;
0078     }
0079
0080     std::vector<const sn*> primitives;
0081     root->collect_prim(primitives);
0082     int  cylinders = 0;
0083     bool outer = false;
0084     bool inner = false;
0085     int  complemented = 0;
0086
0087     for (const sn* primitive : primitives)
0088     {
0089         if (primitive->typecode != CSG_CYLINDER)
0090             continue;
0091
0092         const double* param = primitive->getParam();
0093         const double* aabb = primitive->getAABB();
0094         if (param == nullptr || aabb == nullptr)
0095             continue;
0096
0097         bool phi = close(param[0], startPhi) && close(param[1], deltaPhi);
0098         bool full_phi = close(param[0], 0.) && close(param[1], 0.);
0099         bool bounds = close(aabb[0], -param[3]) && close(aabb[1], -param[3]) && close(aabb[3], param[3]) && close(aabb[4], param[3]);
0100         if (!bounds)
0101             continue;
0102
0103         cylinders++;
0104         complemented += primitive->complement == 1 ? 1 : 0;
0105         outer = outer || (phi && close(param[3], 100.) && primitive->complement == 0);
0106         inner = inner || (full_phi && close(param[3], 50.) && primitive->complement == 1);
0107     }
0108
0109     delete root;
0110
0111     if (cylinders != 2 || !outer || !inner || complemented != 1)
0112     {
0113         std::cerr << "converted tube did not use one wedged outer and one full inner cylinder" << std::endl;
0114         return EXIT_FAILURE;
0115     }
0116
0117     std::cout << "partial-phi G4Tubs converted to a wedged outer and full inner Cylinder" << std::endl;
0118     return EXIT_SUCCESS;
0119 }