eic-opticks/u4/U4SolidMaker.cc

0001 #include <cstring>
0002 #include <csignal>
0003 #include "sstr.h"
0004 #include "ssys.h"
0005 #include "s_bb.h"
0006 #include "NP.hh"
0007
0008 #include "G4SystemOfUnits.hh"
0009 #include "G4Polycone.hh"
0010 #include "G4Sphere.hh"
0011 #include "G4Box.hh"
0012 #include "G4Tubs.hh"
0013 #include "G4Cons.hh"
0014 #include "G4SubtractionSolid.hh"
0015 #include "G4ThreeVector.hh"
0016 #include "G4Orb.hh"
0017 #include "G4UnionSolid.hh"
0018 #include "G4IntersectionSolid.hh"
0019 #include "G4Ellipsoid.hh"
0020 #include "G4MultiUnion.hh"
0021 #include "G4CutTubs.hh"
0022 #include "G4Torus.hh"
0023
0024 using CLHEP::pi ;
0025 using CLHEP::mm ;
0026
0027 #include "U4SolidMaker.hh"
0028 #include "U4SolidTree.hh"
0029 #include "SLOG.hh"
0030
0031 namespace
0032 {
0033     inline const G4VSolid* EnsureVoxelizedMultiUnion(const G4VSolid* solid)
0034     {
0035         const auto* mu = dynamic_cast<const G4MultiUnion*>(solid);
0036         if(mu && mu->GetNumberOfSolids() > 0 && mu->GetVoxels().GetCountOfVoxels() == 0)
0037         {
0038             const_cast<G4MultiUnion*>(mu)->Voxelize();
0039         }
0040         return solid ;
0041     }
0042 }
0043
0044 const plog::Severity U4SolidMaker::LEVEL = SLOG::EnvLevel("U4SolidMaker", "DEBUG");
0045
0046 const char* U4SolidMaker::NAMES = R"LITERAL(
0047 JustOrbOrbUnion
0048 JustOrbOrbIntersection
0049 JustOrbOrbDifference
0050 JustOrb
0051 ThreeOrbUnion
0052 SphereIntersectBox
0053 SphereWithPhiSegment
0054 SphereWithPhiCutDEV
0055 GeneralSphereDEV
0056 SphereWithThetaSegment
0057 AdditionAcrylicConstruction
0058 XJfixtureConstruction
0059 AltXJfixtureConstruction
0060 AltXJfixtureConstructionU
0061 XJanchorConstruction
0062 SJReceiverConstruction
0063 BoxMinusTubs0
0064 BoxMinusTubs1
0065 BoxMinusOrb
0066 UnionOfHemiEllipsoids
0067 PolyconeWithMultipleRmin
0068 PolyconeWithPhiCut
0069 PolyconeWithPhiCutHalf
0070 AnnulusBoxUnion
0071 AnnulusTwoBoxUnion
0072 AnnulusTwoBoxUnionContiguous
0073 AnnulusOtherTwoBoxUnion
0074 AnnulusCrossTwoBoxUnion
0075 AnnulusFourBoxUnion
0076 CylinderFourBoxUnion
0077 BoxFourBoxUnion
0078 BoxCrossTwoBoxUnion
0079 BoxThreeBoxUnion
0080 OrbOrbMultiUnion
0081 OrbGridMultiUnion
0082 BoxGridMultiUnion
0083 BoxFourBoxContiguous
0084 LHCbRichSphMirr
0085 LHCbRichFlatMirr
0086 LocalFastenerAcrylicConstruction
0087 AltLocalFastenerAcrylicConstruction
0088 AnotherLocalFastenerAcrylicConstruction
0089 WaterDistributer
0090 AltWaterDistributer
0091 WaterDistributorPartIIIUnion
0092 OuterReflectorOrbSubtraction
0093 R12860_PMTSolid
0094 LProfileSectorPolycone
0095 )LITERAL";
0096
0097
0098 const char* U4SolidMaker::Name( const char* prefix, unsigned idx ) // static
0099 {
0100     std::stringstream ss ;
0101     ss << prefix << idx ;
0102     std::string s = ss.str();
0103     return strdup(s.c_str()) ;
0104 }
0105
0106
0107 /**
0108 U4SolidMaker::PrimitiveClone
0109 -----------------------------
0110
0111 Create G4VSolid using copy ctor of appropriate EntityType
0112
0113 **/
0114
0115 G4VSolid* U4SolidMaker::PrimitiveClone( const G4VSolid* src, const char* prefix, unsigned idx) // static
0116 {
0117     const char* name = Name(prefix, idx );
0118     G4VSolid* clone = U4SolidTree::PrimitiveClone(src, name);
0119     return clone ;
0120 }
0121
0122 bool U4SolidMaker::StartsWith( const char* n, const char* q ) // static
0123 {
0124     return strlen(q) >= strlen(n) && strncmp(q, n, strlen(n)) == 0 ;
0125 }
0126
0127 bool U4SolidMaker::CanMake(const char* qname) // static
0128 {
0129     bool found = false ;
0130     std::stringstream ss(NAMES) ;
0131     std::string name ;
0132     while (!found && std::getline(ss, name)) if(!name.empty() && StartsWith(name.c_str(), qname)) found = true ;
0133     LOG(LEVEL) << " qname " << qname << " found " << found ;
0134     return found ;
0135 }
0136
0137 const G4VSolid* PolyconeWithMultipleRmin(const char* name);
0138
0139
0140 const G4VSolid* U4SolidMaker::Make(const char* qname)  // static
0141 {
0142     std::string meta ;
0143     return U4SolidMaker::Make(qname, meta);
0144 }
0145 const G4VSolid* U4SolidMaker::Make(const char* qname, std::string& meta )  // static
0146 {
0147     if(strcmp(qname, "NAMES") == 0 )
0148     {
0149         std::cout << NAMES ;
0150         return nullptr ;
0151     }
0152
0153     const G4VSolid* solid = nullptr ;
0154     if(     StartsWith("JustOrbOrbUnion",qname))              solid = JustOrbOrbUnion(qname);
0155     else if(StartsWith("JustOrbOrbIntersection",qname))       solid = JustOrbOrbIntersection(qname);
0156     else if(StartsWith("JustOrbOrbDifference",qname))         solid = JustOrbOrbDifference(qname);
0157     else if(StartsWith("JustOrb",qname))                      solid = JustOrb(qname);
0158     else if(StartsWith("ThreeOrbUnion",qname))                solid = ThreeOrbUnion(qname);
0159     else if(StartsWith("SphereWithPhiSegment",qname))         solid = SphereWithPhiSegment(qname);
0160     else if(StartsWith("SphereWithPhiCutDEV",qname))          solid = SphereWithPhiCutDEV(qname);
0161     else if(StartsWith("GeneralSphereDEV",qname))             solid = GeneralSphereDEV(qname, meta);
0162     else if(StartsWith("SphereWithThetaSegment",qname))       solid = SphereWithThetaSegment(qname);
0163     else if(StartsWith("AdditionAcrylicConstruction",qname))  solid = AdditionAcrylicConstruction(qname);
0164     else if(StartsWith("XJfixtureConstruction", qname))       solid = XJfixtureConstruction(qname);
0165     else if(StartsWith("AltXJfixtureConstructionU", qname))   solid = AltXJfixtureConstructionU(qname);
0166     else if(StartsWith("AltXJfixtureConstruction", qname))    solid = AltXJfixtureConstruction(qname);
0167     else if(StartsWith("XJanchorConstruction", qname))        solid = XJanchorConstruction(qname) ;
0168     else if(StartsWith("SJReceiverConstruction", qname))      solid = SJReceiverConstruction(qname) ;
0169     else if(StartsWith("BoxMinusTubs0",qname))                solid = BoxMinusTubs0(qname);
0170     else if(StartsWith("BoxMinusTubs1",qname))                solid = BoxMinusTubs1(qname);
0171     else if(StartsWith("BoxMinusOrb",qname))                  solid = BoxMinusOrb(qname);
0172     else if(StartsWith("UnionOfHemiEllipsoids", qname))       solid = UnionOfHemiEllipsoids(qname);
0173     else if(StartsWith("PolyconeWithMultipleRmin", qname))    solid = PolyconeWithMultipleRmin(qname) ;
0174     else if(StartsWith("PolyconeWithPhiCutHalf", qname))      solid = PolyconeWithPhiCutHalf(qname) ;
0175     else if(StartsWith("PolyconeWithPhiCut", qname))          solid = PolyconeWithPhiCut(qname) ;
0176     else if(StartsWith("AnnulusBoxUnion", qname))             solid = AnnulusBoxUnion(qname) ;
0177     else if(StartsWith("AnnulusTwoBoxUnion", qname))          solid = AnnulusTwoBoxUnion(qname) ;
0178     else if(StartsWith("AnnulusOtherTwoBoxUnion", qname))     solid = AnnulusOtherTwoBoxUnion(qname) ;
0179     else if(StartsWith("AnnulusCrossTwoBoxUnion", qname))     solid = AnnulusCrossTwoBoxUnion(qname) ;
0180     else if(StartsWith("AnnulusFourBoxUnion", qname))         solid = AnnulusFourBoxUnion(qname) ;
0181     else if(StartsWith("CylinderFourBoxUnion", qname))        solid = CylinderFourBoxUnion(qname) ;
0182     else if(StartsWith("BoxFourBoxUnion", qname))             solid = BoxFourBoxUnion(qname) ;
0183     else if(StartsWith("BoxCrossTwoBoxUnion", qname))         solid = BoxCrossTwoBoxUnion(qname) ;
0184     else if(StartsWith("BoxThreeBoxUnion", qname))            solid = BoxThreeBoxUnion(qname) ;
0185     else if(StartsWith("OrbOrbMultiUnion", qname))            solid = OrbOrbMultiUnion(qname) ;
0186     else if(StartsWith("OrbGridMultiUnion", qname))           solid = OrbGridMultiUnion(qname) ;
0187     else if(StartsWith("BoxGridMultiUnion", qname))           solid = BoxGridMultiUnion(qname) ;
0188     else if(StartsWith("BoxFourBoxContiguous", qname))        solid = BoxFourBoxContiguous(qname) ;
0189     else if(StartsWith("LHCbRichSphMirr", qname))             solid = LHCbRichSphMirr(qname) ;
0190     else if(StartsWith("LHCbRichFlatMirr", qname))            solid = LHCbRichFlatMirr(qname) ;
0191     else if(StartsWith("SphereIntersectBox", qname))          solid = SphereIntersectBox(qname) ;
0192     else if(StartsWith("LocalFastenerAcrylicConstruction", qname)) solid = LocalFastenerAcrylicConstruction(qname) ;
0193     else if(StartsWith("AltLocalFastenerAcrylicConstruction", qname)) solid = AltLocalFastenerAcrylicConstruction(qname) ;
0194     else if(StartsWith("BltLocalFastenerAcrylicConstruction", qname)) solid = BltLocalFastenerAcrylicConstruction(qname) ;
0195     else if(StartsWith("WaterDistributer", qname))                    solid = WaterDistributer(qname) ;
0196     else if(StartsWith("AltWaterDistributer", qname))                 solid = AltWaterDistributer(qname) ;
0197     else if(StartsWith("WaterDistributorPartIIIUnion", qname))        solid = WaterDistributorPartIIIUnion(qname);
0198     else if(StartsWith("OuterReflectorOrbSubtraction", qname))        solid = OuterReflectorOrbSubtraction(qname);
0199     else if(StartsWith("R12860_PMTSolid", qname))                     solid = R12860_PMTSolid(qname);
0200     else if(StartsWith("LProfileSectorPolycone",qname))               solid = LProfileSectorPolycone(qname);
0201     LOG(LEVEL) << " qname " << qname << " solid " << solid ;
0202     LOG_IF(error, solid==nullptr) << " Failed to create solid for qname " << qname << " CHECK U4SolidMaker::Make " ;
0203     return EnsureVoxelizedMultiUnion(solid) ;
0204 }
0205
0206
0207
0208
0209
0210
0211 const G4VSolid* U4SolidMaker::JustOrb(const char* name)  // static
0212 {
0213     return new G4Orb(name, 100.) ;
0214 }
0215
0216
0217 const G4VSolid* U4SolidMaker::JustOrbOrbUnion(       const char* name){ return JustOrbOrb_(name, 'U') ; }
0218 const G4VSolid* U4SolidMaker::JustOrbOrbIntersection(const char* name){ return JustOrbOrb_(name, 'I') ; }
0219 const G4VSolid* U4SolidMaker::JustOrbOrbDifference(  const char* name){ return JustOrbOrb_(name, 'D') ; }
0220
0221 const G4VSolid* U4SolidMaker::JustOrbOrb_(const char* name, char op)  // static
0222 {
0223     G4Orb* l = new G4Orb("l", 100.) ;
0224     G4Orb* r = new G4Orb("r", 100.) ;
0225     G4VSolid* oo  = nullptr ;
0226     switch(op)
0227     {
0228          case 'U': oo = new G4UnionSolid(       name, l, r,  0, G4ThreeVector(0.*mm, 50.*mm, 0.*mm ) ) ; break ;
0229          case 'I': oo = new G4IntersectionSolid(name, l, r,  0, G4ThreeVector(0.*mm, 50.*mm, 0.*mm ) ) ; break ;
0230          case 'D': oo = new G4SubtractionSolid( name, l, r,  0, G4ThreeVector(0.*mm, 50.*mm, 0.*mm ) ) ; break ;
0231     }
0232     return oo ;
0233 }
0234
0235 /**
0236 U4SolidMaker::SphereWithPhiCutDEV
0237 --------------------------------------
0238 **/
0239
0240 const G4VSolid* U4SolidMaker::SphereWithPhiCutDEV(const char* name)  // static
0241 {
0242     double radius    = ssys::getenvfloat("U4SolidMaker_SphereWithPhiCutDEV_radius", 20.f) ;    // mm
0243     double phi_start = ssys::getenvfloat("U4SolidMaker_SphereWithPhiCutDEV_phi_start", 0.f) ;  // units of pi
0244     double phi_delta = ssys::getenvfloat("U4SolidMaker_SphereWithPhiCutDEV_phi_delta", 0.5f) ; // units of pi
0245
0246     G4String pName = name ;
0247     G4double pRmin = 0. ;
0248     G4double pRmax = radius ;
0249     G4double pSPhi = phi_start*pi ;
0250     G4double pDPhi = phi_delta*pi ;
0251     G4double pSTheta = 0. ;
0252     G4double pDTheta = pi ;     // pi: full in theta
0253
0254     return new G4Sphere(pName, pRmin, pRmax, pSPhi, pDPhi, pSTheta, pDTheta );
0255 }
0256
0257
0258
0259
0260 const G4VSolid* U4SolidMaker::GeneralSphereDEV(const char* name, std::string& meta )  // static
0261 {
0262     const char* radiusMode = ssys::getenvvar("U4SolidMaker_GeneralSphereDEV_radiusMode");
0263     double innerRadius = ssys::getenvfloat("U4SolidMaker_GeneralSphereDEV_innerRadius", 50.f) ;    // mm
0264     double outerRadius = ssys::getenvfloat("U4SolidMaker_GeneralSphereDEV_outerRadius", 100.f) ;    // mm
0265
0266     // The two azimuthal angles from the phi cut should be in range 0. to 2. (in units of pi) use XY projection to visualize
0267     const char* phiMode = ssys::getenvvar("U4SolidMaker_GeneralSphereDEV_phiMode");
0268     double phiStart    = ssys::getenvfloat("U4SolidMaker_GeneralSphereDEV_phiStart",    0.f) ;
0269     double phiDelta    = ssys::getenvfloat("U4SolidMaker_GeneralSphereDEV_phiDelta",    2.0f) ;
0270
0271     // The two polar angles from the theta cut should be in range 0. to 1. (in units of pi) use XZ or YZ projections to visualize
0272     const char* thetaMode = ssys::getenvvar("U4SolidMaker_GeneralSphereDEV_thetaMode");
0273     double thetaStart  = ssys::getenvfloat("U4SolidMaker_GeneralSphereDEV_thetaStart",  0.f) ;
0274     double thetaDelta  = ssys::getenvfloat("U4SolidMaker_GeneralSphereDEV_thetaDelta",  0.5f) ;
0275
0276
0277     NP::SetMeta<std::string>( meta, "creator", "U4SolidMaker::GeneralSphereDEV" );
0278     NP::SetMeta<std::string>( meta, "name",    name );
0279
0280     NP::SetMeta<std::string>( meta, "radiusMode", radiusMode );
0281     NP::SetMeta<float>( meta, "innerRadius", float(innerRadius) );
0282     NP::SetMeta<float>( meta, "outerRadius", float(outerRadius) );
0283
0284     NP::SetMeta<std::string>( meta, "phiMode", phiMode );
0285     NP::SetMeta<float>( meta, "phiStart", float(phiStart) );
0286     NP::SetMeta<float>( meta, "phiDelta", float(phiDelta) );
0287
0288     NP::SetMeta<std::string>( meta, "thetaMode", thetaMode );
0289     NP::SetMeta<float>( meta, "thetaStart", float(thetaStart) );
0290     NP::SetMeta<float>( meta, "thetaDelta", float(thetaDelta) );
0291
0292
0293     G4String pName = name ;
0294     G4double pRmin = innerRadius*mm ;
0295     G4double pRmax = outerRadius*mm ;
0296     G4double pSPhi = phiStart*pi ;
0297     G4double pDPhi = phiDelta*pi ;
0298     G4double pSTheta = thetaStart*pi ;
0299     G4double pDTheta = thetaDelta*pi ;
0300
0301     return new G4Sphere(pName, pRmin, pRmax, pSPhi, pDPhi, pSTheta, pDTheta );
0302 }
0303
0304
0305
0306
0307
0308
0309 /**
0310 U4SolidMaker::SphereWithPhiSegment
0311 --------------------------------------
0312
0313 Best way to view phi segment is with XY cross section
0314
0315 phi_start:0 phi_delta:2
0316     full sphere in phi
0317
0318 phi_start:0 phi_delta:0.5
0319     cheese shape : suspect position of the cheese
0320     may differ between Opticks and Geant4
0321
0322
0323 **/
0324
0325 const G4VSolid* U4SolidMaker::SphereWithPhiSegment(const char* name)  // static
0326 {
0327     double phi_start = ssys::getenvfloat("U4SolidMaker_SphereWithPhiSegment_phi_start", 0.f) ;  // units of pi
0328     double phi_delta = ssys::getenvfloat("U4SolidMaker_SphereWithPhiSegment_phi_delta", 0.5f) ; // units of pi
0329
0330     LOG(info)
0331         << " (inputs are scaled by pi) "
0332         << " phi_start  " << phi_start
0333         << " phi_delta  " << phi_delta
0334         << " phi_start+phi_delta  " << phi_start+phi_delta
0335         << " phi_start*pi " << phi_start*pi
0336         << " phi_delta*pi " << phi_delta*pi
0337         << " (phi_start+phi_delta)*pi " << (phi_start+phi_delta)*pi
0338         ;
0339
0340     G4String pName = name ;
0341     G4double pRmin = 0. ;
0342     G4double pRmax = 100. ;
0343     G4double pSPhi = phi_start*pi ;
0344     G4double pDPhi = phi_delta*pi ;
0345     G4double pSTheta = 0. ;
0346     G4double pDTheta = pi ;     // pi: full in theta
0347
0348     return new G4Sphere(pName, pRmin, pRmax, pSPhi, pDPhi, pSTheta, pDTheta );
0349 }
0350
0351
0352
0353
0354
0355
0356 const G4VSolid* U4SolidMaker::SphereWithThetaSegment(const char* name)  // static
0357 {
0358     double theta_start = ssys::getenvfloat("U4SolidMaker_SphereWithThetaSegment_theta_start", 0.f) ;  // units of pi
0359     double theta_delta = ssys::getenvfloat("U4SolidMaker_SphereWithThetaSegment_theta_delta", 0.5f) ; // units of pi
0360
0361     LOG(info)
0362         << " (inputs are scaled by pi) "
0363         << " theta_start  " << theta_start
0364         << " theta_delta  " << theta_delta
0365         << " theta_start+theta_delta  " << theta_start+theta_delta
0366         << " theta_start*pi " << theta_start*pi
0367         << " theta_delta*pi " << theta_delta*pi
0368         << " (theta_start+theta_delta)*pi " << (theta_start+theta_delta)*pi
0369         ;
0370
0371     G4String pName = name ;
0372     G4double pRmin = 0. ;
0373     G4double pRmax = 100. ;
0374     G4double pSPhi = 0.*pi ;
0375     G4double pDPhi = 2.*pi ;
0376     G4double pSTheta = theta_start*pi ;
0377     G4double pDTheta = theta_delta*pi ;  // theta_delta 1. : full in theta
0378
0379     return new G4Sphere(pName, pRmin, pRmax, pSPhi, pDPhi, pSTheta, pDTheta );
0380 }
0381
0382
0383 const G4VSolid* U4SolidMaker::AdditionAcrylicConstruction(const char* name)
0384 {
0385     G4VSolid*   simple             = nullptr ;
0386     G4VSolid*   solidAddition_down = nullptr ;
0387     G4VSolid*   solidAddition_up   = nullptr ;
0388     G4VSolid*   uni_acrylic1       = nullptr ;
0389
0390     double ZNodes3[3];
0391     double RminNodes3[3];
0392     double RmaxNodes3[3];
0393     ZNodes3[0] = 5.7*mm; RminNodes3[0] = 0*mm; RmaxNodes3[0] = 450.*mm;
0394     ZNodes3[1] = 0.0*mm; RminNodes3[1] = 0*mm; RmaxNodes3[1] = 450.*mm;
0395     ZNodes3[2] = -140.0*mm; RminNodes3[2] = 0*mm; RmaxNodes3[2] = 200.*mm;
0396
0397     bool replace_poly = true ;
0398     bool skip_sphere = true ;
0399
0400     simple = new G4Tubs("simple", 0*mm, 450*mm, 200*mm, 0.0*deg,360.0*deg);
0401     solidAddition_down = replace_poly ? simple : new G4Polycone("solidAddition_down",0.0*deg,360.0*deg,3,ZNodes3,RminNodes3,RmaxNodes3);
0402
0403     solidAddition_up = new G4Sphere("solidAddition_up",0*mm,17820*mm,0.0*deg,360.0*deg,0.0*deg,180.*deg);
0404     uni_acrylic1 = skip_sphere ? solidAddition_down : new G4SubtractionSolid("uni_acrylic1",solidAddition_down,solidAddition_up,0,G4ThreeVector(0*mm,0*mm,+17820.0*mm));
0405
0406     G4VSolid*   solidAddition_up1 = nullptr ;
0407     G4VSolid*   solidAddition_up2 = nullptr ;
0408     G4VSolid*   uni_acrylic2      = nullptr ;
0409     G4VSolid*   uni_acrylic3      = nullptr ;
0410
0411     G4VSolid*   uni_acrylic2_initial = nullptr ;
0412
0413     //double zshift = -20*mm ;
0414     double zshift = 0*mm ;
0415
0416     solidAddition_up1 = new G4Tubs("solidAddition_up1",120*mm,208*mm,15.2*mm,0.0*deg,360.0*deg);
0417     uni_acrylic2 = new G4SubtractionSolid("uni_acrylic2",uni_acrylic1,solidAddition_up1,0,G4ThreeVector(0.*mm,0.*mm,zshift));
0418     uni_acrylic2_initial = uni_acrylic2 ;
0419
0420
0421     solidAddition_up2 = new G4Tubs("solidAddition_up2",0,14*mm,52.5*mm,0.0*deg,360.0*deg);
0422
0423     for(int i=0;i<8;i++)
0424     {
0425     uni_acrylic3 = new G4SubtractionSolid("uni_acrylic3",uni_acrylic2,solidAddition_up2,0,G4ThreeVector(164.*cos(i*pi/4)*mm,164.*sin(i*pi/4)*mm,-87.5));
0426     uni_acrylic2 = uni_acrylic3;
0427     }
0428
0429
0430     LOG(info)
0431         << " solidAddition_down " << solidAddition_down
0432         << " solidAddition_up " << solidAddition_up
0433         << " solidAddition_up1 " << solidAddition_up1
0434         << " solidAddition_up2 " << solidAddition_up2
0435         << " uni_acrylic2_initial " << uni_acrylic2_initial
0436         << " uni_acrylic1 " << uni_acrylic1
0437         << " uni_acrylic2 " << uni_acrylic2
0438         << " uni_acrylic3 " << uni_acrylic3
0439         ;
0440
0441     //return solidAddition_down ;  // union of cylinder and cone
0442     //return solidAddition_up1 ;     // pipe cylinder
0443     //return uni_acrylic1 ;
0444     return uni_acrylic2_initial ;
0445 }
0446
0447 /**
0448 U4SolidMaker::ThreeOrbUnion : Checking the cloning of a boolean tree with two levels of transforms
0449 ----------------------------------------------------------------------------------------------------
0450
0451 ThreeOrbUnion0: left-hand form does NOT combine transforms
0452
0453 * this way the transforms are both associated with leaf nodes
0454
0455        abc
0456        /  \
0457       /   (c)  <- transform on leaf
0458      ab
0459     /  \
0460    a   (b)  <- transform on leaf
0461
0462 * all right hand sides are leaf, so no double level of transforms
0463
0464
0465 ThreeOrbUnion1: right-hand form DOES combine two transforms : and trips assert in U4SolidTree::BooleanClone
0466
0467       cab
0468     /    \
0469    c     (ab)   <- transform on operator node
0470           / \
0471          a  (b)  <- transform on leaf
0472
0473 * happens because a boolean combination is on the right hand side of another boolean combination
0474
0475
0476 The righthand edge of 0 is at 120 and 1 is at 130.
0477 They are different because of the combination of the transforms for 1::
0478
0479     XX=120 GEOM=ThreeOrbUnion0_XZ ./xxs.sh
0480     XX=130 GEOM=ThreeOrbUnion1_XZ ./xxs.sh
0481
0482 **/
0483
0484 const G4VSolid* U4SolidMaker::ThreeOrbUnion(const char* name)
0485 {
0486     long mode = sstr::ExtractLong(name, 0);
0487     LOG(LEVEL) << " mode " << mode ;
0488
0489     G4Orb* a = new G4Orb("a", 100. );
0490     G4Orb* b = new G4Orb("b", 100. );
0491     G4Orb* c = new G4Orb("c", 100. );
0492
0493     const G4VSolid* solid = nullptr ;
0494     if( mode == 0 )
0495     {
0496         G4UnionSolid* ab = new G4UnionSolid( "ab", a, b, 0, G4ThreeVector( 10., 0., 0. ) );
0497         G4UnionSolid* ab_c = new G4UnionSolid( "ab_c", ab, c, 0, G4ThreeVector( 20., 0., 0. ) );
0498         solid = ab_c ;
0499     }
0500     else if( mode == 1 )
0501     {
0502         G4UnionSolid* ab = new G4UnionSolid( "ab", a, b, 0, G4ThreeVector( 10., 0., 0. ) );
0503         G4UnionSolid* c_ab = new G4UnionSolid( "c_ab", c, ab, 0, G4ThreeVector( 20., 0., 0. ) );
0504         solid = c_ab ;
0505     }
0506     return solid ;
0507 }
0508
0509 /**
0510 U4SolidMaker::XJfixtureConstruction
0511 -------------------------------------
0512
0513 solidXJfixture             Union
0514    down_uni4               Union
0515       down_uni3            Union
0516          down_uni2         Union
0517              down_uni1     Union
0518                 down1      Tubs
0519                 down2      Box       52.*mm, 0.*mm, 0.*mm
0520              down2         Box      -52.*mm, 0.*mm, 0.*mm
0521          down3             Box       0.*mm,  50.*mm, 0.*mm
0522       down3                Box       0.*mm, -50.*mm, 0.*mm
0523    up_uni                  Union
0524       up1                  Box
0525       up2                  Box       0.*mm, 0.*mm, 13.5*mm     (Up in Z)
0526
0527
0528
0529          up2 is raised by 13.5 to form the thinner in z table top of the altar
0530
0531
0532                 Spurious vertical at 35          In Y box is at 50 +- 15
0533                                                            35    50     65
0534                                                             :
0535                                                             :     :     :                                altar frame              fixture frame
0536                                                             :
0537              -------------+                             +---+---+-+-----+        - - - - - - - - - - - - 18.5+13  =   31.5             6.5     - - - - - -
0538              |            |                             |   :   |      13/2=6.5
0539              +            +                             +   :   + :     :         - - - - - - - - - - -  18.5+6.5 =   25               0.0
0540              |            |                             |   :   |       :
0541              +------------+----------------+-----25-----+---20--+-+-----+         - - - - - - - - - -      8.5+10 =  18.5              -6.5       13+10 = 23
0542              |                                                          5
0543              +    up2                      +                            +       - - - - - - - - - - - - -   8.5+5  = 13.5              -11.5
0544              |                                                          5
0545              +---------+^^^^^^^^^^^^^^^^^^^+^^^^^^^^^^^^^^^^^^+---------+       - - - - - - - - - - - - -             8.5              -16.5    - - - - -
0546                        |                                      |
0547                        |                                     17/2=8.5
0548                        +  up1                                 +                - - - - - - - - - - - - -              0.0              -25.0
0549                        |                                      |
0550                        |                                      |
0551                        +-------------------+-------40---------+            - - - - - - - - - - - - - - - -           -8.5              -33.5
0552
0553                                            |            |    :   |
0554                                            0            25   35  45
0555                                                         |    :   |
0556                                                         |    :   |
0557                                                         |    :   outer tubs
0558                                                         |    :
0559                                                         |    spurious vertical from box edge (why? it is within the tubs ring)
0560                                                         |
0561                                                         inner tubs
0562
0563
0564                Z
0565                |
0566                +-- Y
0567               /
0568              X
0569
0570
0571        Then altar is offset by -25. pushing its top down to 18.5 - 25. = -6.5 in final frame
0572        which is flush with the lower edge of the celtic cross
0573
0574
0575
0576 **/
0577
0578
0579 /**
0580 U4SolidMaker::AltXJfixtureConstruction
0581 ----------------------------------------
0582
0583 Contract this with XJfixtureConstruction : the shape is very nearly the
0584 same but this uses only 3 boxes and 2 tubs rather than 6 boxes and 2 tubs.
0585
0586
0587                         :      65     :
0588                         : 11.5        :
0589         +-----------+---+---+---------+
0590         |           + x - --+ 13      | 23/2    - - - - -
0591         |           +---+---+         |            (23-13)/2 = 10/2
0592         + - - - - - - - | - - - - - - +         - - - - -
0593         |    u          |             | 23/2
0594         |               |             |
0595         +-----+^^^^^^^^^|^^^^^^^^+----+            (23+17)/2  = 40/2
0596               |  l      |        | 17/2
0597               + - - - - | - - - -+              ----
0598               |         |        | 17/2
0599               +---------+---45---+
0600
0601        Z
0602        |
0603        |
0604        +---> Y
0605       /
0606      X
0607                                   ulxoi
0608                       ulxo                    i
0609              ulx              o
0610          ul        x
0611        u    l
0612
0613 **/
0614
0615 const G4VSolid* U4SolidMaker::AltXJfixtureConstruction(const char* name)
0616 {
0617     return AltXJfixtureConstruction_(name, "");
0618 }
0619 const G4VSolid* U4SolidMaker::AltXJfixtureConstructionU(const char* name)
0620 {
0621     return AltXJfixtureConstruction_(name, "ulxo");
0622 }
0623
0624 const G4VSolid* U4SolidMaker::AltXJfixtureConstruction_(const char* name, const char* opt)
0625 {
0626     G4VSolid* u ;
0627     G4VSolid* l ;
0628     G4VSolid* x ;
0629     G4VSolid* o ;
0630     G4VSolid* i ;
0631
0632     G4VSolid* ul ;
0633     G4VSolid* ulx ;
0634     G4VSolid* ulxo ;
0635     G4VSolid* ulxoi ;
0636
0637     // Y is the long (left-right) dimension
0638     G4double l_uncoincide = 1. ;
0639
0640     u = new G4Box("u", 15.*mm, 65.*mm,  23/2.*mm);
0641     l = new G4Box("l", 15.*mm, 40.*mm,  (17+l_uncoincide)/2.*mm);    // increase lbox in half_z by  lbox_uncoincide/2.
0642     ul = new G4UnionSolid("ul", u, l, 0, G4ThreeVector(0.*mm, 0.*mm, (-40.+l_uncoincide)/2*mm ) ) ;
0643
0644     G4double zs = 10/2.*mm ;
0645     x = new G4Box("x", 62.*mm, 11.5*mm, 13/2.*mm);
0646     ulx = new G4UnionSolid("ulx", ul, x, 0, G4ThreeVector(0.*mm, 0.*mm, zs )) ;
0647
0648     o = new G4Tubs("o", 0.*mm, 45.*mm, 13./2*mm, 0.*deg, 360.*deg);
0649     ulxo = new G4UnionSolid("ulxo_CSG_CONTIGUOUS", ulx, o, 0, G4ThreeVector(0.*mm, 0.*mm, zs )) ;
0650     if( strstr(opt, "ulxo"))  return ulxo ;
0651
0652     G4double i_uncoincide = 1. ;
0653     // increase the half_z of subtracted tubs : avoiding upper coincident face
0654     // and raise by the same to keep low edge sub the same : hmm that leaves a thin cross piece from the base of x
0655     // perhaps better to subtract more and get rid of that ?
0656
0657     i = new G4Tubs("i", 0.*mm, 25.*mm, 13./2*mm + i_uncoincide/2.*mm , 0.*deg, 360.*deg);
0658     ulxoi = new G4SubtractionSolid("ulxoi", ulxo, i, 0, G4ThreeVector(0.*mm, 0.*mm, zs+i_uncoincide/2.*mm  )) ;
0659
0660     //return ul ;
0661     //return ulx ;
0662     //return ulxo ;
0663     return ulxoi ;
0664 }
0665 const int U4SolidMaker::XJfixtureConstruction_debug_mode = ssys::getenvint("U4SolidMaker__XJfixtureConstruction_debug_mode", 0 ) ;
0666
0667 /**
0668 U4SolidMaker::XJfixtureConstruction
0669 ------------------------------------
0670
0671
0672 **/
0673
0674 const G4VSolid* U4SolidMaker::XJfixtureConstruction(const char* name)
0675 {
0676     G4VSolid* solidXJfixture_down1;
0677     G4VSolid* solidXJfixture_down2;
0678     G4VSolid* solidXJfixture_down3;
0679     G4VSolid* solidXJfixture_down_uni1;
0680     G4VSolid* solidXJfixture_down_uni2;
0681     G4VSolid* solidXJfixture_down_uni3;
0682     G4VSolid* solidXJfixture_down_uni4;
0683     //G4VSolid* solidXJfixture_down_uni5;
0684
0685     G4VSolid* solidXJfixture_up1;
0686     G4VSolid* solidXJfixture_up2;
0687     G4VSolid* solidXJfixture_up_uni;
0688
0689     G4VSolid* solidXJfixture;
0690
0691 // fixture part
0692     solidXJfixture_down1 = new G4Tubs("solidXJfixture_down1", 25.*mm, 45.*mm, 13./2*mm, 0.*deg, 360.*deg);
0693     solidXJfixture_down2 = new G4Box("solidXJfixture_down2", 10.*mm, 11.5*mm, 13/2.*mm);
0694     solidXJfixture_down_uni1 = new G4UnionSolid("solidXJfixture_down_uni1", solidXJfixture_down1, solidXJfixture_down2, 0, G4ThreeVector(52.*mm, 0.*mm, 0.*mm));
0695     solidXJfixture_down_uni2 = new G4UnionSolid("solidXJfixture_down_uni2", solidXJfixture_down_uni1, solidXJfixture_down2, 0, G4ThreeVector(-52.*mm, 0.*mm, 0.*mm));
0696     solidXJfixture_down3 = new G4Box("solidXJfixture_down3", 15.*mm, 15.*mm, 13/2.*mm);
0697     solidXJfixture_down_uni3 = new G4UnionSolid("solidXJfixture_down_uni3", solidXJfixture_down_uni2, solidXJfixture_down3, 0, G4ThreeVector(0.*mm, 50.*mm, 0.*mm));
0698     solidXJfixture_down_uni4 = new G4UnionSolid("solidXJfixture_down_uni4", solidXJfixture_down_uni3, solidXJfixture_down3, 0, G4ThreeVector(0.*mm, -50.*mm, 0.*mm));
0699
0700     // down_uni4 is celtic-cross shape or uniform z half-thickness 13./2. = 6.5 mm  (shifts in x and y,  not z)
0701
0702 // cover part  : Y is the long dimension
0703     solidXJfixture_up1 = new G4Box("solidXJfixture_up1", 15.*mm, 40.*mm, 17/2.*mm);
0704     solidXJfixture_up2 = new G4Box("solidXJfixture_up2", 15.*mm, 65*mm, 5.*mm);
0705     solidXJfixture_up_uni = new G4UnionSolid("solidXJfixture_up_uni", solidXJfixture_up1, solidXJfixture_up2, 0, G4ThreeVector(0.*mm, 0.*mm, 13.5*mm));
0706
0707     //G4VSolid* new_solidXJfixture_up_uni = Uncoincide_Box_Box_Union( solidXJfixture_up_uni );
0708     //solidXJfixture_up_uni = new_solidXJfixture_up_uni ;
0709
0710     solidXJfixture = new G4UnionSolid("solidXJfixture", solidXJfixture_down_uni4, solidXJfixture_up_uni, 0, G4ThreeVector(0.*mm, 0.*mm, -25.*mm));
0711
0712
0713     // twiddling puts the zero at the altar frame zero
0714     // so would have to offset the placement
0715
0716     G4VSolid* solidXJfixture_twiddle = new G4UnionSolid("solidXJfixture_twiddle", solidXJfixture_up_uni, solidXJfixture_down_uni4, 0, G4ThreeVector(0.*mm, 0.*mm, 25.*mm));
0717     G4VSolid* celtic_cross_sub_altar = new G4SubtractionSolid("solidXJfixture_celtic_cross_sub_altar", solidXJfixture_down_uni4, solidXJfixture_up_uni, 0, G4ThreeVector(0.*mm, 0.*mm, -25.*mm));
0718     G4VSolid* solidXJfixture_split = new G4UnionSolid("solidXJfixture_split", solidXJfixture_down_uni4, solidXJfixture_up_uni, 0, G4ThreeVector(0.*mm, 0.*mm, -50.*mm));
0719
0720
0721     G4VSolid* solid = solidXJfixture ;
0722
0723     int debug_mode = XJfixtureConstruction_debug_mode ;
0724     if( debug_mode > 0 )
0725     {
0726         switch(debug_mode)
0727         {
0728            case  0: solid = solidXJfixture           ; break ;
0729            case  1: solid = solidXJfixture_down1     ; break ;
0730            case  2: solid = solidXJfixture_down2     ; break ;
0731            case  3: solid = solidXJfixture_down_uni1 ; break ;
0732            case  4: solid = solidXJfixture_down_uni2 ; break ;
0733            case  5: solid = solidXJfixture_down3     ; break ;
0734            case  6: solid = solidXJfixture_down_uni3 ; break ;
0735            case  7: solid = solidXJfixture_down_uni4 ; break ;
0736            case  8: solid = solidXJfixture_up1       ; break ;
0737            case  9: solid = solidXJfixture_up2       ; break ;
0738            case 10: solid = solidXJfixture_up_uni    ; break ;
0739            case 11: solid = celtic_cross_sub_altar   ; break ;
0740            case 12: solid = solidXJfixture_split     ; break ;
0741            case 13: solid = solidXJfixture_twiddle   ; break ;
0742         }
0743         LOG(info)
0744             << "U4SolidMaker__XJfixtureConstruction_debug_mode " << debug_mode
0745             << " solid.GetName " << ( solid ? solid->GetName() : "-" )
0746             ;
0747         assert(solid);
0748     }
0749     return solid ;
0750 }
0751
0752 const G4VSolid* U4SolidMaker::AnnulusBoxUnion(const char* name)
0753 {
0754     // do not see spurious intersects
0755     G4VSolid* down1  = new G4Tubs("down1", 25.*mm, 45.*mm, 13./2*mm, 0.*deg, 360.*deg);
0756     G4VSolid* down3 = new G4Box("down3", 15.*mm, 15.*mm, 13/2.*mm);
0757     G4VSolid* uni13 = new G4UnionSolid("uni13", down1, down3, 0, G4ThreeVector(0.*mm, 50.*mm, 0.*mm));
0758     return uni13 ;
0759 }
0760
0761 /**
0762 U4SolidMaker::AnnulusTwoBoxUnion
0763 -----------------------------------
0764
0765 YZ view looks symmetric as both sides of tubs have box-extensions::
0766
0767                                                35       50       65
0768       -65      -45        -25             25    :   45  :        :
0769        +---------+---------+       +       +---------+--:--------+
0770        |         |         |               |         |           |
0771        |         |         |               |         |           |
0772        |         |         |               |         |           |
0773        +---------+---------+               +---------+-----------+
0774
0775     Z
0776     |
0777     +-- Y
0778    /
0779   X
0780
0781
0782 U4SolidMaker::AnnulusTwoBoxUnion
0783
0784
0785                tub_bpy_bny  T
0786                              \
0787      tub_bpy T                bny
0788               \
0789   tub          bpy
0790
0791
0792
0793 For hinting as CSG_CONTIGUOUS or CSG_DISCONTIGUOUS to work
0794 requires setting the inner to zero to avoid the CSG_DIFFERENCE.
0795
0796 **/
0797
0798
0799 const G4VSolid* U4SolidMaker::AnnulusTwoBoxUnion(const char* name)
0800 {
0801     const char* suffix = nullptr ;
0802     if(     strstr(name, "Contiguous"))    suffix = "_CSG_CONTIGUOUS" ;
0803     else if(strstr(name, "Discontiguous")) suffix = "_CSG_DISCONTIGUOUS" ;
0804
0805     bool toplist = strstr(name, "List") != nullptr ;
0806
0807     const char* listname = sstr::Concat("tub_bpy_bny", suffix );
0808
0809     double innerRadius = 0.*mm ;
0810     double bpy_scale_z = 2. ;
0811     double bny_scale_z = 4. ;
0812
0813     LOG(LEVEL)
0814         << " name " << name
0815         << " suffix " << suffix
0816         << " listname " << listname
0817         << " innerRadius " << innerRadius
0818         << " bpy_scale_z " << bpy_scale_z
0819         << " bny_scale_z " << bny_scale_z
0820         ;
0821
0822     G4VSolid* tub  = new G4Tubs("tub", 0.*mm, 45.*mm, 13./2*mm, 0.*deg, 360.*deg);
0823
0824     G4VSolid* bpy = new G4Box("bpy", 15.*mm, 15.*mm, bpy_scale_z*13/2.*mm);
0825     G4VSolid* tub_bpy = new G4UnionSolid(  "tub_bpy", tub, bpy, 0, G4ThreeVector(0.*mm, 50.*mm, 0.*mm));  // +Y
0826
0827     G4VSolid* bny = new G4Box("bny", 15.*mm, 15.*mm, bny_scale_z*13/2.*mm);
0828     G4VSolid* tub_bpy_bny = new G4UnionSolid(listname, tub_bpy, bny, 0, G4ThreeVector(0.*mm, -50.*mm, 0.*mm)); // -Y
0829
0830     double uncoincide = 1. ;
0831     G4VSolid* itu  = new G4Tubs("itu", 0.*mm, 25.*mm, (uncoincide+13./2)*mm, 0.*deg, 360.*deg);
0832     G4VSolid* tub_bpy_bny_itu = new G4SubtractionSolid("tub_bpy_bny_itu", tub_bpy_bny, itu );
0833
0834     return toplist ? tub_bpy_bny : tub_bpy_bny_itu  ;
0835 }
0836
0837
0838
0839 const G4VSolid* U4SolidMaker::AnnulusOtherTwoBoxUnion(const char* name)
0840 {
0841     G4VSolid* down1 = new G4Tubs("down1", 25.*mm, 45.*mm, 13./2*mm, 0.*deg, 360.*deg);
0842     G4VSolid* down2 = new G4Box("down2", 10.*mm, 11.5*mm, 13/2.*mm);
0843     G4VSolid* down_uni1 = new G4UnionSolid("down_uni1", down1    , down2, 0, G4ThreeVector(52.*mm, 0.*mm, 0.*mm));  // +X
0844     G4VSolid* down_uni2 = new G4UnionSolid("down_uni2", down_uni1, down2, 0, G4ThreeVector(-52.*mm, 0.*mm, 0.*mm)); // -X
0845     return down_uni2 ;
0846 }
0847
0848
0849 const G4VSolid* U4SolidMaker::AnnulusCrossTwoBoxUnion(const char* name)
0850 {
0851     G4VSolid* down1 = new G4Tubs("down1", 25.*mm, 45.*mm, 13./2*mm, 0.*deg, 360.*deg);
0852     G4VSolid* down2 = new G4Box("down2", 10.*mm, 11.5*mm, 13/2.*mm);
0853     G4VSolid* down_uni1 = new G4UnionSolid("down_uni1", down1    , down2, 0, G4ThreeVector(52.*mm, 0.*mm, 0.*mm));  // +X
0854
0855     G4VSolid* down3 = new G4Box("down3", 15.*mm, 15.*mm, 13/2.*mm);
0856     G4VSolid* down_uni3 = new G4UnionSolid("down_uni3", down_uni1, down3, 0, G4ThreeVector(0.*mm, 50.*mm, 0.*mm));  // +Y
0857
0858     return down_uni3 ;
0859 }
0860
0861
0862
0863 /**
0864 Not yet managed to see the spurious intersects with a render::
0865
0866     EYE=0,0,2 UP=0,1,0 CAM=1 TMIN=2 ./cxr_geochain.sh
0867
0868 **/
0869
0870 const G4VSolid* U4SolidMaker::AnnulusFourBoxUnion_(const char* name, G4double inner_radius  )
0871 {
0872     // spurious intersects appear in XY cross section but not YZ
0873     G4VSolid* down1 = new G4Tubs("down1", inner_radius, 45.*mm, 13./2*mm, 0.*deg, 360.*deg);
0874     G4VSolid* down2 = new G4Box("down2", 10.*mm, 11.5*mm, 13/2.*mm);
0875     G4VSolid* down_uni1 = new G4UnionSolid("down_uni1", down1    , down2, 0, G4ThreeVector(52.*mm, 0.*mm, 0.*mm));  // +X
0876     G4VSolid* down_uni2 = new G4UnionSolid("down_uni2", down_uni1, down2, 0, G4ThreeVector(-52.*mm, 0.*mm, 0.*mm)); // -X
0877     G4VSolid* down3 = new G4Box("down3", 15.*mm, 15.*mm, 13/2.*mm);
0878     G4VSolid* down_uni3 = new G4UnionSolid("down_uni3", down_uni2, down3, 0, G4ThreeVector(0.*mm, 50.*mm, 0.*mm));  // +Y
0879     G4VSolid* down_uni4 = new G4UnionSolid("down_uni4", down_uni3, down3, 0, G4ThreeVector(0.*mm, -50.*mm, 0.*mm)); // -Y
0880     return down_uni4 ;
0881 }
0882
0883 const G4VSolid* U4SolidMaker::AnnulusFourBoxUnion(const char* name){  return AnnulusFourBoxUnion_(name, 25.*mm );  }
0884 const G4VSolid* U4SolidMaker::CylinderFourBoxUnion(const char* name){ return AnnulusFourBoxUnion_(name,  0.*mm );  }
0885
0886
0887 /**
0888 U4SolidMaker::BoxFourBoxUnion_
0889 ---------------------------------
0890
0891
0892                           +----------+
0893                           | B        |
0894                           |          |
0895                +----------|----------|-------------+
0896                |  A       |          |             |
0897                |          +----------+             |
0898                |                                   |
0899           +---------+                        +------------+
0900           |    |    |                        |     |    C |
0901           |    |    |                        |     |      |
0902           |    |    |                        |     |      |
0903           | D  |    |                        |     |      |
0904           +---------+                        +------------+
0905                |                                   |
0906                |                                   |
0907                |          +----------+             |
0908                |          |          |             |
0909                +----------|----------|-------------+
0910                           |          |
0911                           |       D  |
0912                           +----------+
0913
0914
0915
0916 TODO: switch off balancing and check the impact of pairing order
0917
0918 * eg disjoint unions (like B+C)
0919   although they work on their own they may be implicated with inner boundary spurious
0920
0921
0922
0923 Notice that the tree grows upwards with a new union root
0924 as each prim is union-ed into the combination.
0925
0926                        U
0927                   U       E
0928              U       D
0929         U       C
0930        A  B
0931
0932 :google:`CSG disjoint union`
0933
0934
0935 **/
0936
0937
0938 const G4VSolid* U4SolidMaker::BoxFourBoxUnion_(const char* name, const char* opt )
0939 {
0940     G4VSolid* cbo = new G4Box("cbo", 45.*mm, 45.*mm, 45.*mm );
0941     G4VSolid* xbo = new G4Box("xbo", 10.*mm, 11.5*mm, 13/2.*mm);
0942     G4VSolid* ybo = new G4Box("ybo", 15.*mm, 15.*mm, 13/2.*mm);
0943
0944     bool px = strstr(opt, "+X");
0945     bool nx = strstr(opt, "-X");
0946     bool py = strstr(opt, "+Y");
0947     bool ny = strstr(opt, "-Y");
0948
0949     G4VSolid* comb = nullptr ;
0950     unsigned idx = 0 ;
0951
0952     if(true)
0953     {
0954         comb = PrimitiveClone(cbo,"cbo",idx) ;
0955         idx += 1 ;
0956     }
0957     if(px)
0958     {
0959         comb = new G4UnionSolid(Name("cpx",idx), comb, PrimitiveClone(xbo,"bpx",idx+1), 0, G4ThreeVector(52.*mm, 0.*mm, 0.*mm));  // +X
0960         idx += 2 ;
0961     }
0962     if(nx)
0963     {
0964         comb = new G4UnionSolid(Name("cnx",idx), comb, PrimitiveClone(xbo,"bnx",idx+1), 0, G4ThreeVector(-52.*mm, 0.*mm, 0.*mm)); // -X
0965         idx += 2 ;
0966     }
0967     if(py)
0968     {
0969         comb = new G4UnionSolid(Name("cpy",idx), comb, PrimitiveClone(ybo,"bpy",idx+1), 0, G4ThreeVector(0.*mm, 50.*mm, 0.*mm));  // +Y
0970         idx += 2 ;
0971     }
0972     if(ny)
0973     {
0974         comb = new G4UnionSolid(Name("cny",idx), comb, PrimitiveClone(ybo,"bny",idx+1), 0, G4ThreeVector(0.*mm, -50.*mm, 0.*mm)); // -Y
0975         idx += 2 ;
0976     }
0977
0978     return comb ;
0979 }
0980
0981 const G4VSolid* U4SolidMaker::BoxFourBoxUnion(const char* name ){      return BoxFourBoxUnion_(name, "+X,-X,+Y,-Y") ; }
0982 const G4VSolid* U4SolidMaker::BoxCrossTwoBoxUnion(const char* name ){  return BoxFourBoxUnion_(name, "+X,+Y") ; }
0983 const G4VSolid* U4SolidMaker::BoxThreeBoxUnion(const char* name ){     return BoxFourBoxUnion_(name, "+X,+Y,-Y") ; }
0984
0985
0986
0987 const G4VSolid* U4SolidMaker::BoxFourBoxContiguous_(const char* name, const char* opt )
0988 {
0989     G4VSolid* cbo = new G4Box("cbo", 45.*mm, 45.*mm, 45.*mm );
0990     G4VSolid* xbo = new G4Box("xbo", 10.*mm, 11.5*mm, 13/2.*mm);
0991     G4VSolid* ybo = new G4Box("ybo", 15.*mm, 15.*mm, 13/2.*mm);
0992
0993     bool px = strstr(opt, "+X");
0994     bool nx = strstr(opt, "-X");
0995     bool py = strstr(opt, "+Y");
0996     bool ny = strstr(opt, "-Y");
0997
0998     G4VSolid* item ;
0999     unsigned idx = 0 ;
1000     G4MultiUnion* comb = new G4MultiUnion(name);
1001     G4RotationMatrix rot(0, 0, 0);
1002
1003     if(true)
1004     {
1005         G4ThreeVector pos(0.*mm, 0.*mm, 0.*mm);  // center
1006         G4Transform3D tr(rot, pos);
1007         item = PrimitiveClone(cbo,"cbo",idx) ;
1008         comb->AddNode(*item, tr);
1009         idx+=1 ;
1010     }
1011     if(px)
1012     {
1013         G4ThreeVector pos(52.*mm, 0.*mm, 0.*mm);  // +X
1014         G4Transform3D tr(rot, pos);
1015         item = PrimitiveClone(xbo,"bpx",idx) ;
1016         comb->AddNode(*item, tr);
1017         idx+=1 ;
1018     }
1019     if(nx)
1020     {
1021         G4ThreeVector pos(-52.*mm, 0.*mm, 0.*mm);  // -X
1022         G4Transform3D tr(rot, pos);
1023         item = PrimitiveClone(xbo,"bnx",idx) ;
1024         comb->AddNode(*item, tr);
1025         idx += 1 ;
1026     }
1027     if(py)
1028     {
1029         G4ThreeVector pos(0.*mm, 50.*mm, 0.*mm);  // +Y
1030         G4Transform3D tr(rot, pos);
1031         item = PrimitiveClone(ybo,"bpy",idx) ;
1032         comb->AddNode(*item, tr);
1033         idx += 1 ;
1034     }
1035     if(ny)
1036     {
1037         G4ThreeVector pos(0.*mm, -50.*mm, 0.*mm);  // -Y
1038         G4Transform3D tr(rot, pos);
1039         item = PrimitiveClone(ybo,"bny",idx) ;
1040         comb->AddNode(*item, tr);
1041         idx += 1 ;
1042     }
1043     return comb ;
1044 }
1045
1046
1047 const G4VSolid* U4SolidMaker::BoxFourBoxContiguous(const char* name ){     return BoxFourBoxContiguous_(name, "-X,+X,+Y,-Y") ; }
1048
1049
1050
1051
1052
1053
1054 /**
1055 U4SolidMaker::Uncoincide_Box_Box_Union
1056 ----------------------------------------
1057
1058 To avoid coincidence need to expand the smaller box into the larger
1059 without changing the position the lower edge.
1060 Hence increase the half-size in Z from *hz* to *new_hz* and
1061 simultaneously shift upwards by the same amount (*zoff*)
1062 to keep the lower edge at same z position::
1063
1064
1065        +hz + uncoincide - - - - - - - +~~~~~~~~~+ - -    zoff + new_hz  - - -
1066                                       |         |
1067                                       |         |
1068        +hz  +---------+ - - - - - - - | - - - - | - - - - - - - - - - - - - -
1069             |         |               |         |
1070             |         |               |         |
1071             |         |               |         |
1072             |         |               |_________|        zoff
1073             |         |               |         |
1074         0 --|---------| - - - - - - - - - - - - - - - - - - - - - - - - - - -
1075             |         |               |         |
1076             |         |               |         |
1077             |         |               |         |
1078             |         |               |         |
1079             |         |               |         |
1080        -hz  +---------+ - - - - - - - +---------+ - - -  zoff - new_hz  - - -
1081
1082
1083
1084 Line equations::
1085
1086       hz + uncoincide = zoff + new_hz
1087
1088       -hz             = zoff - new_hz
1089
1090 Add them::
1091
1092      zoff = uncoincide/2
1093
1094 Subtract them::
1095
1096      new_hz = hz + uncoincide/2
1097
1098
1099
1100
1101
1102 up2 is raised by 13.5 to form the thinner in z table top of the altar
1103
1104
1105          +---------------------------+     5mm                     - - -  8.5 + 10 = 18.5
1106          |         up2               |  - - - -   13.5  = 8.5+5
1107          +-----+---------------+-----+
1108                |               |   17/2 = 8.5mm
1109                |   up1         |   - - - -
1110                |               |
1111                +---------------+
1112                                           10 mm thin top of altar,
1113                                           17 mm thicker bottom of altar
1114
1115
1116
1117
1118 HMM: the sign of the change to the translation depends on
1119 whether the smaller_box (which needs to grow into the larger) is on
1120 the rhs of the combination which has the translation applied to it
1121
1122
1123 **/
1124
1125 G4VSolid* U4SolidMaker::Uncoincide_Box_Box_Union( const G4VSolid* bbu  )  // static
1126 {
1127     LOG(info) << " bbu.GetName " << bbu->GetName() ;
1128
1129     const G4BooleanSolid* bs = dynamic_cast<const G4BooleanSolid*>(bbu) ;
1130
1131     const G4VSolid* a = bs->GetConstituentSolid(0) ;
1132     const G4VSolid* _b = bs->GetConstituentSolid(1) ;
1133     const G4DisplacedSolid* _b_disp = dynamic_cast<const G4DisplacedSolid*>(_b) ;
1134     G4ThreeVector b_tla = _b_disp->GetObjectTranslation();
1135     const G4VSolid* b = _b_disp->GetConstituentMovedSolid()  ;
1136
1137     LOG(info) << " a.GetName " << a->GetName() ;
1138     LOG(info) << " _b.GetName " << _b->GetName() ;
1139     LOG(info) << " _b_disp.GetName " << _b_disp->GetName() ;
1140     LOG(info) << " b.GetName " << b->GetName() ;
1141     LOG(info) << " b_tla " << Desc(&b_tla) ;
1142
1143     const G4Box* a_box = dynamic_cast<const G4Box*>(a);
1144     const G4Box* b_box = dynamic_cast<const G4Box*>(b);
1145     LOG(info) << " a_box " << Desc(a_box) ;
1146     LOG(info) << " b_box " << Desc(b_box) ;
1147
1148     G4ThreeVector new_tla(b_tla);
1149
1150     std::string new_name = bs->GetName()  ;
1151
1152     G4Box* new_a = new G4Box(a->GetName(), a_box->GetXHalfLength(), a_box->GetYHalfLength(), a_box->GetZHalfLength() );
1153     G4Box* new_b = new G4Box(b->GetName(), b_box->GetXHalfLength(), b_box->GetYHalfLength(), b_box->GetZHalfLength() );
1154
1155     int shift_axis = OneAxis(&b_tla);
1156     LOG(info) << " shift_axis " << shift_axis ;
1157
1158
1159
1160     enum { A, B, UNKNOWN } ;
1161     int smaller = UNKNOWN ;
1162
1163     for(int axis=0 ; axis < 3 ; axis++)
1164     {
1165         if(axis == shift_axis) continue ;   //
1166         double ah = HalfLength(a_box, axis);
1167         double bh = HalfLength(b_box, axis);
1168         if(ah == bh) continue ;    // ignore equal axes
1169         smaller = ah < bh ? A : B ;
1170     }
1171
1172     LOG(info) << " smaller " << smaller ;
1173
1174
1175     double uncoincide = 1.*mm ;
1176     if(smaller != UNKNOWN )
1177     {
1178         G4Box* smaller_box = smaller == A ? new_a : new_b ;
1179         LOG(info) << " smaller_box.GetName " << smaller_box->GetName() ;
1180
1181         ChangeBoxHalfLength( smaller_box, shift_axis, uncoincide/2. );
1182
1183         ChangeThreeVector(   &new_tla ,   shift_axis, uncoincide/2. ) ;
1184     }
1185     else
1186     {
1187         LOG(fatal) << " failed to uncoincide " ;
1188     }
1189
1190     G4UnionSolid* new_bbu = new G4UnionSolid( new_name, new_a, new_b, 0, new_tla );
1191     return new_bbu ;
1192 }
1193
1194
1195 double U4SolidMaker::HalfLength( const G4Box* box, int axis )
1196 {
1197     double value = 0. ;
1198     switch(axis)
1199     {
1200         case X: value = box->GetXHalfLength() ; break ;
1201         case Y: value = box->GetYHalfLength() ; break ;
1202         case Z: value = box->GetZHalfLength() ; break ;
1203     }
1204     return value ;
1205 }
1206
1207 void U4SolidMaker::ChangeThreeVector( G4ThreeVector* v, int axis, double delta )
1208 {
1209     if( v == nullptr ) return ;
1210     switch(axis)
1211     {
1212         case X: v->setX(v->x() + delta)  ; break ;
1213         case Y: v->setY(v->y() + delta)  ; break ;
1214         case Z: v->setZ(v->z() + delta)  ; break ;
1215     }
1216 }
1217 void U4SolidMaker::ChangeBoxHalfLength( G4Box* box, int axis, double delta )
1218 {
1219     switch(axis)
1220     {
1221         case X: box->SetXHalfLength(box->GetXHalfLength() + delta)  ; break ;
1222         case Y: box->SetYHalfLength(box->GetYHalfLength() + delta)  ; break ;
1223         case Z: box->SetZHalfLength(box->GetZHalfLength() + delta)  ; break ;
1224     }
1225 }
1226
1227
1228
1229
1230
1231
1232 int U4SolidMaker::OneAxis( const G4ThreeVector* v )
1233 {
1234     double x = v ? v->x() : 0. ;
1235     double y = v ? v->y() : 0. ;
1236     double z = v ? v->z() : 0. ;
1237     int axis = ERR ;
1238     if( x != 0. && y == 0. && z == 0. ) axis = X ;
1239     if( x == 0. && y != 0. && z == 0. ) axis = Y ;
1240     if( x == 0. && y == 0. && z != 0. ) axis = Z ;
1241     return axis ;
1242 }
1243
1244
1245
1246 std::string U4SolidMaker::Desc( const G4Box* box )
1247 {
1248     std::stringstream ss ;
1249     ss
1250        << "("
1251        << std::fixed << std::setw(10) << std::setprecision(3) << box->GetXHalfLength() << " "
1252        << std::fixed << std::setw(10) << std::setprecision(3) << box->GetYHalfLength() << " "
1253        << std::fixed << std::setw(10) << std::setprecision(3) << box->GetZHalfLength()
1254        << ") "
1255        << box->GetName()
1256        ;
1257     std::string s = ss.str();
1258     return s ;
1259 }
1260
1261 std::string U4SolidMaker::Desc( const G4ThreeVector* v )
1262 {
1263     std::stringstream ss ;
1264     ss
1265        << "("
1266        << std::fixed << std::setw(10) << std::setprecision(3) << (v ? v->x() : 0. ) << " "
1267        << std::fixed << std::setw(10) << std::setprecision(3) << (v ? v->y() : 0. )  << " "
1268        << std::fixed << std::setw(10) << std::setprecision(3) << (v ? v->z() : 0. )
1269        << ") "
1270        ;
1271     std::string s = ss.str();
1272     return s ;
1273 }
1274
1275
1276
1277
1278 /**
1279 U4SolidMaker::XJanchorConstruction
1280 -----------------------------------
1281
1282 Observed spurious Geant4 intersects on the line between the Tubs and the Cons::
1283
1284     solidXJanchor          G4UnionSolid
1285
1286         sub                G4SubtractionSolid      (subtract big sphere from cone)
1287               down         G4Cons
1288               ball         G4Sphere
1289
1290         up                 G4Tubs
1291
1292
1293
1294
1295         +-------------------------+--------------------------+         - - - - - -
1296          \                        .                         /                             10.0
1297           .                       +                        .           - - - - - -                   - - - - - -
1298            \                      .                       /                               10.0            |
1299             +---------+^^^^^^^^^^^.^^^^^^^^^^^^+---------+             - - - - - -                       16.5
1300                       |           .            |                                    13/2 = 6.5            |
1301                       +           .            +                       - - - - - -                   - - - - - -
1302                       |           .            |                                    13/2 = 6.5
1303                       +-----------.------------+         |   |         - - - - - -
1304                                   0           25        47  73
1305
1306
1307      FIX :
1308          increase tubs hz by uncoincide/2
1309          shift upwards uncoincide/2 (-> low edge stays same) by shifting down less
1310
1311 **/
1312
1313 const int U4SolidMaker::XJanchorConstruction_debug_mode = ssys::getenvint("U4SolidMaker__XJanchorConstruction_debug_mode", 0 ) ;
1314
1315 const G4VSolid* U4SolidMaker::XJanchorConstruction(const char* name)
1316 {
1317     bool do_uncoincide = false ;
1318     bool do_noball = false ;
1319
1320     switch(XJanchorConstruction_debug_mode)
1321     {
1322         case 0: do_uncoincide = false ;  do_noball = false ; break ;
1323         case 1: do_uncoincide = true  ;  do_noball = false ; break ;
1324         case 2: do_uncoincide = false ;  do_noball = true  ; break ;
1325         case 3: do_uncoincide = true  ;  do_noball = true  ; break ;
1326     }
1327
1328     double uncoincide = do_uncoincide ? 1. : 0. ;
1329
1330     LOG(info)
1331         << " U4SolidMaker__XJanchorConstruction_debug_mode " << XJanchorConstruction_debug_mode
1332         << " do_uncoincide " << do_uncoincide
1333         << " uncoincide " << uncoincide
1334         << " do_noball " << do_noball
1335         ;
1336
1337
1338     G4VSolid* solidXJanchor_up;
1339     G4VSolid* solidXJanchor_down;
1340     G4VSolid* solidXJanchor_ball;
1341
1342     solidXJanchor_up   = new G4Tubs("solidXJanchor_up", 0.*mm, 25.*mm, (13.+uncoincide)/2*mm, 0.*deg, 360.*deg);
1343     solidXJanchor_down = new G4Cons("solidXJanchor_down", 0.*mm, 47.*mm, 0.*mm, 73.*mm, 10.*mm, 0.*deg, 360.*deg);   // to subtract the ball
1344     solidXJanchor_ball = new G4Sphere("solidXJanchor_ball", 0.*mm, 17820.*mm, 0.*deg, 360.*deg, 0.*deg, 180.*deg);
1345
1346     G4SubtractionSolid* solidXJanchor_sub_ = new G4SubtractionSolid("solidXJanchor_sub",solidXJanchor_down, solidXJanchor_ball, 0, G4ThreeVector(0.*mm, 0*mm,  17820.*mm));
1347     G4VSolid* solidXJanchor_sub = do_noball ? solidXJanchor_down : (G4VSolid*)solidXJanchor_sub_ ;
1348
1349     G4UnionSolid* solidXJanchor = new G4UnionSolid("solidXJanchor",solidXJanchor_sub, solidXJanchor_up, 0, G4ThreeVector(0.*mm, 0*mm,(-16.5 + uncoincide/2)*mm));
1350
1351     return solidXJanchor ;
1352 }
1353
1354
1355
1356
1357 const G4VSolid* U4SolidMaker::SJReceiverConstruction(const char* name)
1358 {
1359 #ifdef DIRTY
1360     G4VSolid* solidSJReceiver_up ;
1361     G4VSolid* solidSJReceiver_down ;
1362     G4VSolid* solidSJReceiver_box ;
1363     G4VSolid* solidSJReceiver_ball ;
1364
1365     //solidSJReceiver_up   = new G4Cons("solidSJReceiver_up", 0.*mm, 31.7*mm, 0.*mm, 25*mm, 13./2*mm, 0.*deg, 360.0*deg);
1366     solidSJReceiver_up = new G4Tubs("solidXJanchor_up", 0.*mm, 25.*mm, 13./2*mm, 0.*deg, 360.*deg);
1367     solidSJReceiver_down = new G4Cons("solidSJReceiver_down", 0.*mm, 73.*mm, 0.*mm, 47.*mm, 10.*mm, 0.*deg, 360.*deg);   // to subtract the ball
1368     //solidSJReceiver_down = new G4Cons("solidSJReceiver_down", 0.*mm, 47.*mm, 0.*mm, 60.*mm, 5.*mm, 0.*deg, 360.*deg); // original size
1369     solidSJReceiver_box = new G4Box("solidSJReceiver_box", 17780.*mm, 17780.*mm, 17780.*mm);
1370     solidSJReceiver_ball = new G4Sphere("solidSJReceiver_ball", 0.*mm, 17700.*mm, 0.*deg, 360.*deg, 0.*deg, 180.*deg);
1371     G4SubtractionSolid* solidSphere_sub = new G4SubtractionSolid("solidSphere_sub", solidSJReceiver_box, solidSJReceiver_ball);
1372     G4SubtractionSolid* solidSJReceiver_sub = new G4SubtractionSolid("solidSJReceiver_sub",solidSJReceiver_down, solidSJReceiver_ball, 0, G4ThreeVector(0.*mm, 0*mm,  17699.938*mm));
1373     G4UnionSolid* solidSJReceiver = new G4UnionSolid("solidSJReceiver",solidSJReceiver_sub, solidSJReceiver_up, 0, G4ThreeVector(0.*mm, 0*mm,16.5*mm));
1374 #else
1375     G4VSolid* solidSJReceiver_up ;
1376     G4VSolid* solidSJReceiver_down ;
1377     G4VSolid* solidSJReceiver_ball ;
1378
1379     solidSJReceiver_up = new G4Tubs("solidXJanchor_up", 0.*mm, 25.*mm, 13./2*mm, 0.*deg, 360.*deg);
1380     solidSJReceiver_down = new G4Cons("solidSJReceiver_down", 0.*mm, 73.*mm, 0.*mm, 47.*mm, 10.*mm, 0.*deg, 360.*deg);   // to subtract the ball
1381     solidSJReceiver_ball = new G4Sphere("solidSJReceiver_ball", 0.*mm, 17700.*mm, 0.*deg, 360.*deg, 0.*deg, 180.*deg);
1382     G4SubtractionSolid* solidSJReceiver_sub = new G4SubtractionSolid("solidSJReceiver_sub",solidSJReceiver_down, solidSJReceiver_ball, 0, G4ThreeVector(0.*mm, 0*mm,  17699.938*mm));
1383     G4UnionSolid* solidSJReceiver = new G4UnionSolid("solidSJReceiver",solidSJReceiver_sub, solidSJReceiver_up, 0, G4ThreeVector(0.*mm, 0*mm,16.5*mm));
1384 #endif
1385     return solidSJReceiver ;
1386 }
1387
1388
1389
1390
1391 const G4VSolid* U4SolidMaker::BoxMinusTubs0(const char* name)  // is afflicted
1392 {
1393     double tubs_hz = 15.2*mm ;
1394     double zshift = 0*mm ;
1395     G4VSolid* box   = new G4Box("box",  250*mm, 250*mm, 100*mm );
1396     G4VSolid* tubs =  new G4Tubs("tubs",120*mm,208*mm,tubs_hz,0.0*deg,360.0*deg);
1397     G4VSolid* box_minus_tubs = new G4SubtractionSolid(name,box,tubs,0,G4ThreeVector(0.*mm,0.*mm,zshift));
1398     return box_minus_tubs ;
1399 }
1400
1401 const G4VSolid* U4SolidMaker::BoxMinusTubs1(const char* name)
1402 {
1403     double tubs_hz = 15.2*mm ;
1404     G4VSolid* box   = new G4Box("box",  250*mm, 250*mm, 100*mm );
1405     G4VSolid* tubs =  new G4Tubs("tubs",120*mm,208*mm,tubs_hz,0.0*deg,360.0*deg);
1406     G4VSolid* box_minus_tubs = new G4SubtractionSolid(name,box,tubs);
1407     return box_minus_tubs ;
1408 }
1409
1410 const G4VSolid* U4SolidMaker::BoxMinusOrb(const char* name)
1411 {
1412     double radius = ssys::getenvfloat("U4SolidMaker_BoxMinusOrb_radius", 130.f) ;
1413
1414     double sx     = ssys::getenvfloat("U4SolidMaker_BoxMinusOrb_sx", 100.f) ;
1415     double sy     = ssys::getenvfloat("U4SolidMaker_BoxMinusOrb_sy", 100.f) ;
1416     double sz     = ssys::getenvfloat("U4SolidMaker_BoxMinusOrb_sz", 100.f) ;
1417
1418     double dx     = ssys::getenvfloat("U4SolidMaker_BoxMinusOrb_dx", 0.f) ;
1419     double dy     = ssys::getenvfloat("U4SolidMaker_BoxMinusOrb_dy", 0.f) ;
1420     double dz     = ssys::getenvfloat("U4SolidMaker_BoxMinusOrb_dz", 0.f) ;
1421
1422     G4VSolid* box = new G4Box("box",  sx, sy, sz );
1423     G4VSolid* orb = new G4Orb("orb",  radius );
1424
1425     G4VSolid* box_minus_orb = new G4SubtractionSolid(name,box,orb,nullptr, G4ThreeVector(dx, dy, dz) );
1426     return box_minus_orb ;
1427 }
1428
1429
1430 const G4VSolid* U4SolidMaker::PolyconeWithMultipleRmin(const char* name)
1431 {
1432     double ZUpper4[4];
1433     double RminUpper4[4];
1434     double RmaxUpper4[4];
1435     ZUpper4[0] = 0*mm;   RminUpper4[0] = 43.*mm; RmaxUpper4[0] = 195.*mm;
1436     ZUpper4[1] = -15*mm; RminUpper4[1] = 43.*mm; RmaxUpper4[1] = 195.*mm;
1437     ZUpper4[2] = -15*mm; RminUpper4[2] = 55.5*mm; RmaxUpper4[2] = 70.*mm;
1438     ZUpper4[3] = -101*mm; RminUpper4[3] = 55.5*mm; RmaxUpper4[3] = 70.*mm;
1439
1440     G4VSolid* base_steel = new G4Polycone("base_steel",0.0*deg,360.0*deg,4,ZUpper4,RminUpper4,RmaxUpper4);
1441     return base_steel ;
1442 }
1443
1444 /**
1445 PolyconeWithPhiCut
1446 ------------------
1447
1448 phiStart=30 phiDelta=30 selects wedge in middle of +X+Y quadrant::
1449
1450    SOLID=PolyconeWithPhiCut EYE=0,0,1000  UP=0,1,0 ~/o/u4/tests/U4SolidMakerTest.sh
1451
1452
1453           Y 30. + +
1454           |  . 30+ .
1455           | .+  .
1456           |.+.   30
1457           +-------X
1458          /
1459         /
1460        Z
1461
1462 **/
1463
1464 const G4VSolid* U4SolidMaker::PolyconeWithPhiCut(const char* name)
1465 {
1466     const int N = 2 ;
1467
1468     double ZPlane[N] = {  -50*mm ,  50*mm } ;
1469     double Rmin[N]   = {    0*mm ,   0*mm } ;
1470     double Rmax[N]   = {  200*mm , 200*mm } ;
1471
1472     double phiStart = 30.0*deg ;
1473     double phiDelta = 30.0*deg ;
1474     G4VSolid* polycone = new G4Polycone(name,phiStart,phiDelta,Z,ZPlane,Rmin,Rmax);
1475     return polycone ;
1476 }
1477
1478 /**
1479 U4SolidMaker::PolyconeWithPhiCutHalf
1480 --------------------------------------
1481
1482 Confirmed that 0->pi selects +Y hemi with::
1483
1484     SOLID=PolyconeWithPhiCutHalf EYE=0,0,1000  UP=0,1,0 ~/o/u4/tests/U4SolidMakerTest.sh
1485
1486 Convert this with::
1487
1488     GEOM # set GEOM to LocalPolyconeWithPhiCutHalf
1489     ~/o/g4cx/tests/G4CX_U4TreeCreateCSGFoundryTest.sh
1490
1491 Viz with::
1492
1493      cxr_min.sh
1494         # from axial viewpoint get full circle of cylinder,
1495         # from less axial get changing geom with a crevasse forming
1496         # in the middle and getting larger as change angle until
1497         # get expected half circle
1498
1499         # Cannot pin down any halfspace bug, so suspect its coming from
1500         # the CSG unbounded handling
1501         #
1502         # WIP: try to reproduce issue on CPU with CSG/tests/csg_intersect_prim_test.sh
1503
1504 **/
1505
1506 const G4VSolid* U4SolidMaker::PolyconeWithPhiCutHalf(const char* name)
1507 {
1508     const int N = 2 ;
1509     double ZPlane[N] = {  -50*mm ,  50*mm } ;
1510     double Rmin[N]   = {    0*mm ,   0*mm } ;
1511     double Rmax[N]   = {  200*mm , 200*mm } ;
1512
1513     double phiStart = 0.0*deg ;
1514     double phiDelta = 180.0*deg ;
1515     G4VSolid* polycone = new G4Polycone(name,phiStart,phiDelta,Z,ZPlane,Rmin,Rmax);
1516     return polycone ;
1517 }
1518
1519
1520
1521
1522
1523
1524
1525
1526 const G4VSolid* U4SolidMaker::UnionOfHemiEllipsoids(const char* name )
1527 {
1528     assert( strstr( name, "UnionOfHemiEllipsoids" ) != nullptr );
1529
1530     std::vector<long> vals ;
1531     sstr::Extract(vals, name);
1532     long iz = vals.size() > 0 ? vals[0] : 0 ;
1533
1534     std::cout
1535         << "U4SolidMaker::UnionOfHemiEllipsoids"
1536         << " name " << name
1537         << " vals.size " << vals.size()
1538         << " iz " << iz
1539         << std::endl
1540         ;
1541
1542
1543     double rx = 150. ;
1544     double ry = 150. ;
1545     double rz = 100. ;
1546
1547     double z2 = rz ;
1548     double z1 = 0. ;
1549     double z0 = -rz ;
1550
1551     G4VSolid* upper = new G4Ellipsoid("upper", rx, ry, rz, z1, z2 );
1552     G4VSolid* lower = new G4Ellipsoid("lower", rx, ry, rz, z0, z1 );
1553
1554     G4VSolid* solid = nullptr ;
1555     if( iz == 0 )
1556     {
1557         solid = new G4UnionSolid(name, upper, lower );
1558     }
1559     else
1560     {
1561         double zoffset = double(iz) ;
1562         G4ThreeVector tlate(0., 0., zoffset) ;
1563         solid = new G4UnionSolid(name, upper, lower, nullptr, tlate );
1564     }
1565     return solid ;
1566 }
1567
1568
1569 /**
1570 U4SolidMaker::OrbOrbMultiUnion
1571 -------------------------------
1572
1573 OrbOrbMultiUnion0
1574     single Orb within MultiUnion
1575 OrbOrbMultiUnion1
1576     three Orb at (-100,0,100) x positions within MultiUnion
1577 OrbOrbMultiUnion2
1578     five Orb at (-200,-100,0,100,200) x positions within MultiUnion
1579 OrbOrbMultiUnion{N}
1580     2N+1 Orb at {-N, -N+1, ..., 0, 1, .... N}*100 x positions within MultiUnion
1581
1582 **/
1583
1584 const G4VSolid* U4SolidMaker::OrbOrbMultiUnion(const char* name )
1585 {
1586     long num = sstr::ExtractLong(name, 1) ;
1587     G4MultiUnion* comb = new G4MultiUnion(name);
1588
1589     G4RotationMatrix rot(0, 0, 0);
1590     G4Orb* orb = new G4Orb("Orb", 50.) ;
1591
1592     for(long i=-num ; i <= num ; i++)
1593     {
1594        G4ThreeVector tlate( double(i)*100.*mm, 0.*mm, 0.*mm);
1595        G4Transform3D tr(rot, tlate) ;
1596        comb->AddNode( *orb, tr );
1597     }
1598
1599     return comb ;
1600 }
1601
1602
1603
1604
1605
1606 /**
1607 U4SolidMaker::GridMultiUnion_
1608 ---------------------------------
1609
1610 (nx, ny, nz) of (3,3,3) yields a 7x7x7 grid from eg::
1611
1612     nx : -3, -2, -1, 0, +1, +2, +3
1613
1614 **/
1615
1616
1617 const G4VSolid* U4SolidMaker::GridMultiUnion_(const char* name, G4VSolid* item, double gridspace, int nx, int ny, int nz )
1618 {
1619     G4MultiUnion* grid = new G4MultiUnion(name);
1620
1621     for(int i=-nx ; i <= nx ; i++ )
1622     for(int j=-ny ; j <= ny ; j++ )
1623     for(int k=-nz ; k <= nz ; k++ )
1624     {
1625         G4ThreeVector pos(double(i)*gridspace*mm, double(j)*gridspace*mm, double(k)*gridspace*mm );
1626         LOG(info) << pos ;
1627
1628         G4RotationMatrix rot(0, 0, 0);
1629         G4Transform3D tr(rot, pos);
1630         grid->AddNode(*item, tr);
1631
1632     }
1633     // The Geant4 header in G4MultiUnion.hh explicitly says Voxelize() must be called once before navigation use.
1634     grid->Voxelize();
1635     return grid ;
1636 }
1637
1638 /**
1639 U4SolidMaker::OrbGridMultiUnion
1640 ----------------------------------
1641
1642                     :       :
1643                     :       :
1644           +---+   +-:-+   +-:-+
1645           |   |   | : |   | : |
1646           +---+   +-:-+   +-:-+
1647                     :       :
1648           +---+   +-:-+   +-:-+
1649           |   |   | 0 |   | : |
1650           +---+   +-:-+   +-:-+
1651                     :       :
1652           +---+   +-:-+   +-:-+
1653           |   |   | : |   | : |
1654           +---+   +-:-+   +-:-+
1655                     :       :
1656                   : : :     :
1657             -radius:+radius
1658                     :
1659
1660
1661 * when there is no overlap  radius < gridscale , only the middle Orb gets intersects
1662
1663 **/
1664
1665 const G4VSolid* U4SolidMaker::OrbGridMultiUnion(const char* name)
1666 {
1667     std::vector<long> vals ;
1668     sstr::Extract(vals, name);
1669
1670     assert( vals.size() == 2 );
1671
1672     double radius(vals[0]) ;
1673     double gridscale(vals[1]) ;
1674
1675     LOG(info)
1676         << " name " << name
1677         << " radius " << radius
1678         << " gridscale " << gridscale
1679         ;
1680
1681     G4VSolid* item = new G4Orb("orb", radius*mm );
1682
1683     int nx = 3 ;
1684     int ny = 3 ;
1685     int nz = 3 ;
1686
1687     return GridMultiUnion_(name, item, gridscale, nx, ny, nz );
1688 }
1689
1690
1691 /**
1692 U4SolidMaker::BoxGridMultiUnion
1693 --------------------------------
1694
1695 eg BoxGridMultiUnion10:30_YX
1696
1697 halfside   10 mm
1698 gridscale  30
1699
1700 **/
1701
1702
1703 const G4VSolid* U4SolidMaker::BoxGridMultiUnion(const char* name)
1704 {
1705     std::vector<long> vals ;
1706     sstr::Extract(vals, name);
1707
1708     assert( vals.size() == 2 );
1709
1710     double halfside(vals[0]) ;
1711     double gridscale(vals[1]) ;
1712
1713     LOG(info)
1714         << " name " << name
1715         << " halfside " << halfside
1716         << " gridscale " << gridscale
1717         ;
1718
1719     G4VSolid* item = new G4Box("box", halfside*mm, halfside*mm, halfside*mm );
1720
1721     int nx = 3 ;
1722     int ny = 3 ;
1723     int nz = 3 ;
1724
1725     return GridMultiUnion_(name, item, gridscale, nx, ny, nz );
1726 }
1727
1728
1729
1730
1731
1732
1733
1734
1735 /**
1736 U4SolidMaker::LHCbRichSphMirr
1737 --------------------------------
1738
1739
1740                                Z
1741                                |
1742
1743               |+--------------|+
1744                |               |
1745                |               |
1746                |               |
1747                |               |
1748                |               |  1500/2
1749                |               |
1750                |               |
1751                |               |
1752                + - - - + - - - +
1753                |               |
1754                |               |
1755                |               |
1756                |               |  1500/2
1757                |               |
1758                |               |
1759                |               |
1760                |               |
1761               |+--------------|+
1762                :               :
1763              3650.0           3651.5
1764
1765
1766                    ------->  X
1767
1768
1769
1770
1771
1772 **/
1773
1774
1775 const G4VSolid* U4SolidMaker::LHCbRichSphMirr(const char* qname)  // static
1776 {
1777     long mode = sstr::ExtractLong(qname, 0);
1778     LOG(LEVEL) << " qname " << qname << " mode " << mode ;
1779
1780    // /Users/blyth/liyu/Rich_Simplified/include/SphMirrorGeometryParameters.hh
1781     //Parameters for simplified rich1 sph mirror.
1782     const G4double InnerRadius = 3650.0 *  CLHEP::mm;
1783     const G4double Thickness   = 1.5 * CLHEP::mm;
1784     const G4double OuterRadius =  InnerRadius + Thickness ;
1785
1786     const G4double SegmentSizeX = 1500.0 * CLHEP::mm;  // SCB : misnammed should be Z not X
1787     const G4double SegmentSizeZ = SegmentSizeX  ;
1788     const G4double SegmentSizeY = 650.0  * CLHEP::mm;
1789
1790
1791     // To account for spherical geometry with straight edges, the angluar sizes are increased by a factor and then boolen subtraction made.
1792     const G4double DeltaExtendedFactor = 1.2;
1793     const G4double DeltaTheta =  2 * asin(SegmentSizeX/(2.0*InnerRadius)) * DeltaExtendedFactor *  CLHEP::rad;
1794     const G4double DeltaPhi =  2 * asin(SegmentSizeY/(2.0*InnerRadius) )* DeltaExtendedFactor * CLHEP::rad;
1795     const G4double ThetaStart = (0.5 * CLHEP::pi * CLHEP::rad)  - (0.5 * DeltaTheta) ;
1796     const G4double PhiStart   =  -0.5 * DeltaPhi ;
1797     const G4double SubLargeSize = 5000 *  CLHEP::mm;
1798     const G4double SubHalfSizeX = SubLargeSize;
1799     const G4double SubHalfSizeY = SubLargeSize;
1800     const G4double SubHalfSizeZ = SubLargeSize;
1801     const G4double SubPosX = 0.0 *  CLHEP::mm;
1802     const G4double SubPosY = 0.0 *  CLHEP::mm;
1803     const G4double SubUpPosY = (0.5 * SegmentSizeY) + SubHalfSizeY;
1804     const G4double SubDownPosY = -1.0 * SubUpPosY ;
1805     const G4double SubPosZ = 0.0 * CLHEP::mm;
1806     const G4double SubRightPosZ =(0.5 * SegmentSizeX) +  SubHalfSizeZ;
1807     const G4double SubLeftPosZ = -1.0* SubRightPosZ;
1808
1809
1810
1811    // /Users/blyth/liyu/Rich_Simplified/src/RichTbLHCbR1or.cc
1812
1813     G4Sphere* SphFullSph = new G4Sphere ("SphFullSphDEV",InnerRadius,
1814                        OuterRadius,PhiStart,
1815                                            DeltaPhi, ThetaStart,
1816                                            DeltaTheta);
1817     G4Box * SphSubBox = new G4Box("SphSubBox",SubHalfSizeX,
1818                     SubHalfSizeY, SubHalfSizeZ);
1819     G4RotationMatrix SubBoxRot;
1820     G4ThreeVector  SubBoxTPos ( SubPosX,SubUpPosY,SubPosZ);
1821     G4ThreeVector  SubBoxBPos ( SubPosX,SubDownPosY,SubPosZ);
1822     G4ThreeVector  SubBoxLPos ( SubPosX,SubPosY,SubLeftPosZ);
1823     G4ThreeVector  SubBoxRPos ( SubPosX,SubPosY,SubRightPosZ);
1824     G4Transform3D  SubBoxTTrans(  SubBoxRot, SubBoxTPos);
1825     G4Transform3D  SubBoxBTrans(  SubBoxRot, SubBoxBPos);
1826     G4Transform3D  SubBoxLTrans(  SubBoxRot, SubBoxLPos);
1827     G4Transform3D  SubBoxRTrans(  SubBoxRot, SubBoxRPos);
1828     G4SubtractionSolid* TSph = new G4SubtractionSolid("SphTSph",SphFullSph,
1829                                             SphSubBox , SubBoxTTrans);
1830    G4SubtractionSolid* BSph = new G4SubtractionSolid("SphBSphBox",TSph,
1831                                             SphSubBox , SubBoxBTrans);
1832    G4SubtractionSolid* LSph = new G4SubtractionSolid("SphLSphBox",BSph,
1833                                             SphSubBox , SubBoxLTrans);
1834    G4SubtractionSolid* RSph = new G4SubtractionSolid("SphRSphBox",LSph,
1835                                             SphSubBox , SubBoxRTrans);
1836
1837
1838     const G4VSolid* simple = nullptr ;
1839     if( mode == 1 )
1840     {
1841         G4double MiddleRadius = (InnerRadius + OuterRadius)/2. ;
1842
1843         G4Orb* inner = new G4Orb("inner", InnerRadius );
1844         G4Orb* outer = new G4Orb("outer", OuterRadius );
1845         G4SubtractionSolid* shell = new G4SubtractionSolid("shell", outer, inner );
1846
1847         G4double sag_max = SagittaMax( InnerRadius, OuterRadius, SegmentSizeY, SegmentSizeZ );
1848         G4double SagMax = 80.*mm ;
1849         bool sag_max_expect = SagMax > sag_max ;
1850         assert( sag_max_expect );
1851         if(!sag_max_expect) std::raise(SIGINT);
1852         G4double FullDepthX = Thickness + SagMax ;
1853
1854         G4Box* cutbox = new G4Box("cutbox", FullDepthX/2. , SegmentSizeY/2. , SegmentSizeZ/2. );
1855         simple = new G4IntersectionSolid("simple_CSG_EXBB", shell, cutbox, 0, G4ThreeVector(MiddleRadius, 0., 0.) ) ;
1856     }
1857
1858
1859     const G4VSolid* solid = nullptr ;
1860     switch(mode)
1861     {
1862        case 0: solid = RSph   ; break ;
1863        case 1: solid = simple ; break ;
1864     }
1865     assert(solid);
1866     return solid  ;
1867 }
1868 /**
1869 U4SolidMaker::Sagitta
1870 ---------------------
1871
1872
1873
1874                   .
1875             .     .        .
1876         .         s           .
1877       /           .            \
1878    A +------------+-------------+ B
1879           l/2     |C    l/2
1880        .          |           .
1881                  r-s
1882           .       |      .
1883         r         |         r
1884               .   |   .
1885                   |
1886                   +
1887                   O
1888
1889 s: sagitta or depth of the arc
1890
1891 l: full length of chord across base of arc
1892
1893 r: radius
1894
1895
1896 Pythagoras for triangle AOC::
1897
1898      (r-s)^2 + (l/2)^2 = r^2
1899
1900        r - s = sqrt( r^2 - (l/2)^2 )
1901
1902            s = r - sqrt( r^2 - (l/2)^2 )
1903
1904 **/
1905
1906
1907 G4double U4SolidMaker::Sagitta( G4double radius, G4double chord )
1908 {
1909     G4double sagitta = radius - sqrt( radius*radius - chord*chord/4. ) ;
1910     return sagitta ;
1911 }
1912
1913 G4double U4SolidMaker::SagittaMax( G4double radius, G4double sy, G4double sz )
1914 {
1915     G4double sag_sy = Sagitta( radius, sy );
1916     G4double sag_sz = Sagitta( radius, sz );
1917     G4double sag_max = std::max( sag_sy, sag_sz );
1918     return sag_max ;
1919 }
1920
1921 G4double U4SolidMaker::SagittaMax( G4double InnerRadius, G4double OuterRadius,  G4double sy, G4double sz )
1922 {
1923     G4double sag_inner = SagittaMax( InnerRadius, sy, sz );
1924     G4double sag_outer = SagittaMax( OuterRadius, sy, sz );
1925     G4double sag_max = std::max( sag_inner, sag_outer );
1926
1927     LOG(LEVEL)
1928         << " InnerRadius " << InnerRadius
1929         << " OuterRadius " << OuterRadius
1930         << " sag_max " << sag_max
1931         << std::endl
1932         << " sag_inner : SaggitaMax(" << InnerRadius << "," <<  sy << "," << sz << ")" << sag_inner
1933         << std::endl
1934         << " sag_outer : SaggitaMax(" << OuterRadius << "," <<  sy << "," << sz << ")" << sag_outer
1935         << std::endl
1936         ;
1937
1938     return sag_max ;
1939 }
1940
1941
1942 /**
1943 U4SolidMaker::LHCbRichFlatMirr
1944 ---------------------------------
1945
1946
1947                   Z
1948                   |
1949                   |    Y
1950                   |   /
1951                   |  /
1952                   | /
1953                   |/                                                                                               |
1954                   +-------> X - - - - - - - - - - - - 5,000,000 mm  - - - - - - - - - - - - - - - - - - - - - - - -+
1955                  .                                                                                                /
1956                 .
1957                .
1958              -Y
1959
1960
1961       Rectangle is formed close to X axis, deltaPhi (in XY plane) is symmetric around phi=0.
1962
1963
1964
1965
1966                                Z
1967                                |
1968
1969               |+--------------|+
1970                |               |
1971                |               |
1972                |               |
1973                |               |
1974                |               |  1480/2
1975                |               |
1976                |               |
1977                |               |
1978                + - - - + - - - +
1979                |               |
1980                |               |
1981                |               |
1982                |               |  1480/2
1983                |               |
1984                |               |
1985                |               |
1986                |               |
1987               |+--------------|+
1988
1989             5,000,000     5,000,006
1990
1991                    ------->  X
1992
1993
1994    sagitta
1995
1996          l^2        ( 1480/2 )^2
1997          ---   =     -------------   = 0.01369
1998          8r            8*5,000,000
1999
2000
2001    intersect box
2002        center    :  ( 5000003, 0 , 0 )
2003        halfsides :  (       4, 880/2,   1480/2 )
2004
2005
2006                     4 in X is more than enough to contain the sagitta
2007
2008 **/
2009
2010 const G4VSolid* U4SolidMaker::LHCbRichFlatMirr(const char* qname)  // static
2011 {
2012     long mode = sstr::ExtractLong(qname, 0);
2013     LOG(LEVEL) << " qname " << qname << " mode " << mode ;
2014
2015
2016     // /Users/blyth/Rich_Simplified/Rich_Simplified/include/RichTbR1FlatMirrorGeometryParameters.hh
2017
2018     //  :.,+30s/RichTbR1FlatMirr//g
2019
2020     //Parameters for simplified rich1 flat mirror.
2021     const G4double InnerRadius = 5000.0 *  CLHEP::m;
2022     const G4double Thickness   = 6.0 * CLHEP::mm;
2023     const G4double OuterRadius =  InnerRadius + Thickness ;
2024     const G4double MiddleRadius = (InnerRadius + OuterRadius)/2. ;
2025
2026
2027     const G4double SegmentSizeX = 1480.0 * CLHEP::mm;  // <-- should be Z (not X)
2028     const G4double SegmentSizeZ = SegmentSizeX  ;
2029
2030     const G4double SegmentSizeY = 880.0 * CLHEP::mm;
2031
2032    /**
2033    From xxs.sh looking at XZ, XY, YZ planes
2034
2035          X range :  5M -> 5M + 6 mm
2036          Y range : -440 -> 440 mm
2037          Z range : -740 -> 740 mm
2038    **/
2039
2040
2041    /*
2042     const G4double BotInLHCbPosY  = 337.90 * CLHEP::mm;
2043     const G4double BotInLHCbPosZ  = 1323.31 * CLHEP::mm;
2044     const G4double VertTilt = 0.25656 * CLHEP:: rad;
2045     const G4double RotX = -1.0 * VertTilt;
2046     const G4double RotY = (0.5 * CLHEP::pi * CLHEP::rad);
2047    */
2048
2049     const G4double DeltaTheta = asin(SegmentSizeX/InnerRadius) * CLHEP::rad;  // NO NEED FOR THE asin AS EXTREMELY SMALL ANGLES
2050     const G4double DeltaPhi   = asin(SegmentSizeY/InnerRadius) * CLHEP::rad;
2051     const G4double StartTheta = (0.5 * CLHEP::pi * CLHEP::rad)  - (0.5*  DeltaTheta);
2052     const G4double StartPhi   = -0.5 * DeltaPhi;
2053
2054     LOG(LEVEL)
2055         << " DeltaTheta " << DeltaTheta
2056         << " SegmentSizeX " << SegmentSizeX
2057         << " InnerRadius " << InnerRadius
2058         << " SegmentSizeX/InnerRadius " << SegmentSizeX/InnerRadius
2059         << " asin(SegmentSizeX/InnerRadius) " << asin(SegmentSizeX/InnerRadius)
2060         << " CLHEP::rad " << CLHEP::rad
2061         << " StartTheta " << StartTheta
2062         ;
2063
2064
2065     LOG(LEVEL)
2066         << " DeltaPhi " << DeltaPhi
2067         << " SegmentSizeY " << SegmentSizeY
2068         << " InnerRadius " << InnerRadius
2069         << " SegmentSizeY/InnerRadius " << SegmentSizeY/InnerRadius
2070         << " asin(SegmentSizeY/InnerRadius) " << asin(SegmentSizeY/InnerRadius)
2071         << " CLHEP::rad " << CLHEP::rad
2072         << " StartPhi " << StartPhi
2073         ;
2074
2075
2076     /*
2077     const G4double InLHCbPosX = 0.0* CLHEP::mm;
2078     const G4double InLHCbPosY = BotInLHCbPosY +
2079                                                  0.5 * SegmentSizeY * cos (VertTilt) +
2080                                                  OuterRadius * sin ( VertTilt );
2081      */
2082
2083
2084     // /Users/blyth/Rich_Simplified/Rich_Simplified/src/RichTbLHCbR1FlatMirror.cc
2085
2086     G4Sphere* RichTbR1FlatFull = new G4Sphere ("RichTbR1FlatFullDEV",InnerRadius,OuterRadius,StartPhi,DeltaPhi, StartTheta,DeltaTheta);
2087
2088
2089
2090     G4VSolid* solid = nullptr ;
2091     if( mode == 0 )
2092     {
2093         solid = RichTbR1FlatFull ;
2094     }
2095     else
2096     {
2097         LOG(LEVEL) << " constructing the flat mirror using intersection box rather than attemping to describe with miniscule thetacut and phicut angles"  ;
2098
2099
2100         G4double sag_max = SagittaMax( InnerRadius, OuterRadius, SegmentSizeY, SegmentSizeZ );
2101         G4double SagMax = 1.*mm ;
2102         assert( SagMax > sag_max );
2103         G4double FullDepthX = Thickness + SagMax ;
2104
2105
2106         LOG(LEVEL)
2107             << " sag_max " << sag_max
2108             << " SagMax " << SagMax
2109             << " FullDepthX  " << FullDepthX
2110             ;
2111
2112         G4Orb* outer = new G4Orb("outer",OuterRadius);
2113         G4Orb* inner = new G4Orb("inner",InnerRadius);
2114         G4SubtractionSolid* shell = new G4SubtractionSolid("shell", outer, inner );
2115         G4Box* box = new G4Box("box", FullDepthX/2. , SegmentSizeY/2. , SegmentSizeZ/2. );
2116
2117         std::string rootname = qname ;
2118         rootname += "_CSG_EXBB" ;
2119         solid = new G4IntersectionSolid( rootname, shell, box, 0, G4ThreeVector( MiddleRadius, 0., 0. ));
2120     }
2121
2122     return solid ;
2123 }
2124
2125
2126
2127
2128 const G4VSolid* U4SolidMaker::SphereIntersectBox(const char* qname)  // static
2129 {
2130     G4double phiStart = 0. ;
2131     G4double phiDelta = 2. ;
2132     G4double thetaStart = 0. ;
2133     G4double thetaDelta = 1.  ;
2134     G4Sphere* sph = new G4Sphere("sph", 95.*mm, 105.*mm,  phiStart*CLHEP::pi, phiDelta*CLHEP::pi, thetaStart*CLHEP::pi, thetaDelta*CLHEP::pi );
2135     G4Box* box = new G4Box("box", 20.*mm, 20.*mm, 20.*mm );
2136     G4IntersectionSolid* sph_box = new G4IntersectionSolid("sph_box_CSG_EXBB", sph, box, 0, G4ThreeVector(0.*mm, 0.*mm, 100.*mm ) );
2137     return sph_box ;
2138 }
2139
2140 /**
2141 U4SolidMaker::LocalFastenerAcrylicConstruction
2142 ------------------------------------------------
2143
2144 The tree is grown upwards::
2145
2146
2147                                 uni1
2148                               /     \
2149                           uni1       screw
2150                         /      \
2151                     uni1        screw
2152                   /     \
2153                uni1      screw
2154               /    \
2155        IonRing     screw
2156
2157
2158
2159 **/
2160
2161
2162
2163 const G4VSolid* U4SolidMaker::LocalFastenerAcrylicConstruction(const char* name) // static
2164 {
2165     const char* PREFIX = "LocalFastenerAcrylicConstruction" ;
2166     assert( sstr::StartsWith(name,PREFIX ));
2167     int num_column = strlen(name) > strlen(PREFIX) ? std::atoi( name + strlen(PREFIX) ) : 8 ;
2168
2169     LOG(info)
2170         << " name " <<  ( name ? name : "-" )
2171         << " num_column " << num_column
2172         ;
2173
2174
2175     const char* screw_name = "screw_HINT_LISTNODE_PRIM_CONTIGUOUS" ;
2176     //const char* screw_name = "screw_HINT_LISTNODE_PRIM_DISCONTIGUOUS" ;
2177
2178     G4VSolid* uni_Addition(nullptr);
2179     {
2180         G4Tubs *IonRing = new G4Tubs("IonRing",123*mm,206.2*mm,7*mm,0.0*deg,360.0*deg);
2181         G4Tubs* screw = new G4Tubs(screw_name,0,13*mm,50.*mm,0.0*deg,360.0*deg);
2182         uni_Addition = IonRing;
2183         for(int i=0;i<num_column;i++)
2184         {
2185             G4ThreeVector tlate(164.*cos(i*pi/4)*mm, 164.*sin(i*pi/4)*mm,-65.0*mm);
2186             G4UnionSolid* uni1 = new G4UnionSolid("uni1",uni_Addition, screw, 0, tlate);
2187             uni_Addition = uni1;
2188         }
2189     }
2190     return uni_Addition ;
2191 }
2192
2193
2194 /**
2195 U4SolidMaker::AltLocalFastenerAcrylicConstruction
2196 --------------------------------------------------
2197
2198          uni1
2199         /     \
2200      IonRing    muni
2201
2202
2203 **/
2204
2205 const G4VSolid* U4SolidMaker::AltLocalFastenerAcrylicConstruction(const char* name) // static
2206 {
2207     [[maybe_unused]] const char* PREFIX = "AltLocalFastenerAcrylicConstruction" ;
2208     assert( sstr::StartsWith(name,PREFIX ));
2209     long num_column = sstr::ExtractLong(name, 1) ;
2210
2211     LOG(info)
2212         << " name " <<  ( name ? name : "-" )
2213         << " num_column " << num_column
2214         ;
2215
2216     assert( num_column > 0 );
2217
2218     G4Tubs* IonRing = new G4Tubs("IonRing",123*mm,206.2*mm,7*mm,0.0*deg,360.0*deg);
2219
2220     G4MultiUnion* muni = new G4MultiUnion(name);
2221     G4Tubs* screw = new G4Tubs("screw",0,13*mm,50.*mm,0.0*deg,360.0*deg);
2222
2223     G4RotationMatrix rot(0, 0, 0);
2224     for(long i=0;i<num_column;i++)
2225     {
2226        G4ThreeVector tlate(164.*cos(i*pi/4)*mm, 164.*sin(i*pi/4)*mm,-65.0*mm);
2227        G4Transform3D tr(rot, tlate) ;
2228        muni->AddNode( *screw, tr );
2229     }
2230
2231     G4UnionSolid* uni1 = new G4UnionSolid(name, IonRing, muni, 0, G4ThreeVector(0.,0.,0.));
2232     return uni1 ;
2233 }
2234
2235
2236 /**
2237 U4SolidMaker::BltLocalFastenerAcrylicConstruction
2238 --------------------------------------------------
2239
2240 Was trying to form a boolean tree  with the nodes destined
2241 to become constituents of the listnode within a separate G4UnionSolid tree.
2242 But it is not so easy to split off the screws like that, as they all need translation
2243 including the 0th.
2244
2245 The reason for trying to do that is because the inner radius on the IonRing tubs
2246 means have to keep that separate.
2247
2248 **/
2249
2250 const G4VSolid* U4SolidMaker::BltLocalFastenerAcrylicConstruction(const char* name) // static
2251 {
2252     [[maybe_unused]] const char* PREFIX = "BltLocalFastenerAcrylicConstruction" ;
2253     assert( sstr::StartsWith(name,PREFIX ));
2254     const size_t num_column = sstr::ExtractSize(name, 1) ;
2255
2256     LOG(info)
2257         << " name " <<  ( name ? name : "-" )
2258         << " num_column " << num_column
2259         ;
2260
2261     assert( num_column > 0 );
2262
2263     G4Tubs* IonRing = new G4Tubs("IonRing",123*mm,206.2*mm,7*mm,0.0*deg,360.0*deg);
2264
2265     G4Tubs* screw = new G4Tubs("screw",0,13*mm,50.*mm,0.0*deg,360.0*deg);
2266
2267     std::vector<G4ThreeVector> tlate(num_column, G4ThreeVector(0, 0, 0));
2268     for(size_t i=0;i<num_column;i++) tlate[i] = G4ThreeVector(164.*cos(i*pi/4)*mm, 164.*sin(i*pi/4)*mm,-65.0*mm);
2269
2270     G4VSolid* muni = screw ;
2271     for(size_t i=1 ; i < num_column ; i++) muni = new G4UnionSolid( name, muni, screw, 0, tlate[i] ) ;
2272
2273     G4UnionSolid* uni1 = new G4UnionSolid(name, IonRing, muni, 0, tlate[0] );
2274     return uni1 ;
2275 }
2276
2277
2278
2279 const G4VSolid* U4SolidMaker::WaterDistributer(const char* name_) // static
2280 {
2281     G4MultiUnion* multiUnion = new G4MultiUnion(name_);
2282
2283     double TubeR_type_I = 38 * mm;
2284     double TubeR_type_II = 70 * mm;
2285     double holeR = TubeR_type_I ;
2286
2287
2288
2289     G4double distance = 120 *cm;
2290     G4double TubeR_I = TubeR_type_I;
2291     G4double TubeR_II = TubeR_type_II;
2292     G4double TubeR_III = TubeR_type_I;
2293
2294     G4double holeRadius = holeR;
2295
2296     G4String name = multiUnion->GetName() ;
2297     G4String name_I = name + "_I" ;
2298     G4String name_II = name + "_II" ;
2299     G4String name_III = name + "_III" ;
2300
2301     G4VSolid* solid_I = WaterDistributerHelper(name_I.c_str(),distance, TubeR_I);
2302     G4VSolid* solid_II = WaterDistributerHelper(name_II.c_str(),distance, TubeR_II);
2303     G4VSolid* solid_III = WaterDistributerHelper(name_III.c_str(),distance, TubeR_III);
2304
2305     G4ThreeVector positionI(0, 23.5*cm, 0);
2306     G4ThreeVector positionII(0, 0, 0);
2307     G4ThreeVector positionIII(0, - 23.5*cm, 0);
2308
2309     G4RotationMatrix* rotation = new G4RotationMatrix();
2310
2311     G4Transform3D trI = G4Transform3D(*rotation, positionI);
2312     G4Transform3D trII = G4Transform3D(*rotation, positionII);
2313     G4Transform3D trIII = G4Transform3D(*rotation, positionIII);
2314
2315     multiUnion->AddNode(*solid_II, trII);
2316     multiUnion->AddNode(*solid_I, trI);
2317     multiUnion->AddNode(*solid_III, trIII);
2318
2319     for (int i = 0; i < 8; i++) {
2320         // Create circular hole geometry
2321         G4Tubs* hole = new G4Tubs("Hole_" + std::to_string(i), 0, holeRadius, (23.5) * cm, 0, 360 * deg);
2322
2323         G4double theta = i * 45 * deg;
2324         G4RotationMatrix* rotation = new G4RotationMatrix();
2325         rotation->rotateX(90 * deg);
2326
2327         G4ThreeVector position(-distance * cos(theta), 0, distance * sin(theta));
2328         G4Transform3D tr = G4Transform3D(*rotation, position);
2329
2330         multiUnion->AddNode(*hole, tr);
2331     }
2332
2333     multiUnion->Voxelize();
2334
2335     return multiUnion ;
2336 }
2337
2338
2339
2340
2341
2342
2343
2344
2345
2346
2347
2348
2349
2350 /**
2351 U4SolidMaker::WaterDistributerHelper
2352 -------------------------------------
2353
2354
2355 ::
2356
2357                _________
2358               / _______ \
2359              / /       \ \
2360             / /         \ \
2361             | |         | |
2362             | |         | |
2363             | |         | |
2364             \ \         / /
2365              \ \_______/ /
2366               \_________/
2367
2368 **/
2369
2370
2371
2372
2373 G4VSolid* U4SolidMaker::WaterDistributerHelper(const char* name, double distance, double TubeR) // static
2374 {
2375     G4double rMax = TubeR;
2376     G4double rMin = 0;
2377     G4double dz = distance * tan(22.5 * deg);
2378     G4double sPhi = 0.0;
2379     G4double dPhi = 2 * M_PI;
2380     //G4double holeRadius = 33 *mm;
2381
2382     G4MultiUnion* multiUnion = new G4MultiUnion(name);
2383
2384     // Rotating and combining eight cylinders
2385     G4double rotationAngle = 22.5 * deg;
2386     for (int i = 0; i < 8; i++) {
2387         // Calculate the normal vector of the cutting plane
2388         G4double theta_i = rotationAngle;
2389         G4double theta_j = rotationAngle;
2390         G4ThreeVector pLowNorm(sin(theta_i), 0, -cos(theta_j));
2391         G4ThreeVector pHighNorm(sin(theta_i), 0, cos(theta_j));
2392
2393         G4CutTubs* cutTube = new G4CutTubs(
2394             "CutTube_" + std::to_string(i), rMin, rMax, dz, sPhi, dPhi, pLowNorm, pHighNorm);
2395
2396         //G4Tubs* hole = new G4Tubs("Hole_" + std::to_string(i), 0, holeRadius, TubeR + 1 * cm, 0, 360 * deg);
2397
2398         G4double theta = i * 45 * deg;
2399         G4RotationMatrix* rotation = new G4RotationMatrix();
2400         rotation->rotateY(45 * i * deg);
2401
2402         G4RotationMatrix* rotation_hole = new G4RotationMatrix();
2403         rotation_hole->rotateX(90 * deg);
2404
2405         G4ThreeVector position(-distance * cos(theta), 0, distance * sin(theta));
2406         G4Transform3D tr = G4Transform3D(*rotation, position);
2407
2408
2409         //G4SubtractionSolid* cutTubeWithHole = new G4SubtractionSolid(
2410         //    "CutTubeWithHole_" + std::to_string(i), cutTube, hole, rotation_hole, G4ThreeVector(0, 0, 0));
2411
2412
2413         multiUnion->AddNode(*cutTube, tr);
2414     }
2415
2416     // Complete the construction of MultiUnion (must call!).
2417     multiUnion->Voxelize();
2418
2419     return multiUnion ;
2420 }
2421
2422
2423
2424
2425 /**
2426 U4SolidMaker::AltWaterDistributer
2427 ----------------------------------
2428
2429 One multiunion of everything works and is simpler
2430 but polygonization taking forever so unusable.
2431
2432 **/
2433
2434 const G4VSolid* U4SolidMaker::AltWaterDistributer(const char* name_) // static
2435 {
2436     G4MultiUnion* multiUnion = new G4MultiUnion(name_);
2437
2438     double TubeR_type_I = 38 * mm;
2439     double TubeR_type_II = 70 * mm;
2440     double holeR = TubeR_type_I ;
2441
2442     G4double distance = 120 *cm;
2443     G4double TubeR_I = TubeR_type_I;
2444     G4double TubeR_II = TubeR_type_II;
2445     G4double TubeR_III = TubeR_type_I;
2446
2447     G4double holeRadius = holeR;
2448
2449     G4double y_positionI = 23.5*cm ;
2450     G4double y_positionII = 0 ;
2451     G4double y_positionIII = -23.5*cm ;
2452
2453     bool do_I = true ;
2454     bool do_II = false ;
2455     bool do_III = false ;
2456
2457     if(do_I)   AltWaterDistributerHelper(multiUnion, distance, TubeR_I,   y_positionI   );
2458     if(do_II)  AltWaterDistributerHelper(multiUnion, distance, TubeR_II,  y_positionII  );
2459     if(do_III) AltWaterDistributerHelper(multiUnion, distance, TubeR_III, y_positionIII );
2460
2461     for (int i = 0; i < 8; i++) {
2462         // Create circular hole geometry : MISNAMED THEY ARE CROSS PIECES IN Y-DIRECTION
2463         G4Tubs* hole = new G4Tubs("Hole_" + std::to_string(i), 0, holeRadius, (23.5) * cm, 0, 360 * deg);
2464
2465         G4double theta = i * 45 * deg;
2466         G4RotationMatrix* rotation = new G4RotationMatrix();
2467         rotation->rotateX(90 * deg);  // G4Tubs starts on Z-axis, rotate around X by 90 puts on Y axis
2468
2469         G4ThreeVector position(-distance * cos(theta), 0, distance * sin(theta));
2470         G4Transform3D tr = G4Transform3D(*rotation, position);
2471
2472         multiUnion->AddNode(*hole, tr);
2473     }
2474
2475     multiUnion->Voxelize();
2476
2477     return multiUnion ;
2478 }
2479
2480
2481 void U4SolidMaker::AltWaterDistributerHelper(G4MultiUnion* multiUnion, double distance, double TubeR, double yoffset ) // static
2482 {
2483     G4double rMax = TubeR;
2484     G4double rMin = 0;
2485     G4double dz = distance * tan(22.5 * deg);
2486     G4double sPhi = 0.0;
2487     G4double dPhi = 2 * M_PI;
2488     //G4double holeRadius = 33 *mm;
2489
2490     // Rotating and combining eight cylinders
2491     G4double rotationAngle = 22.5 * deg;
2492     for (int i = 0; i < 8; i++) {
2493         // Calculate the normal vector of the cutting plane
2494         G4double theta_i = rotationAngle;
2495         G4double theta_j = rotationAngle;
2496         G4ThreeVector pLowNorm(sin(theta_i), 0, -cos(theta_j));
2497         G4ThreeVector pHighNorm(sin(theta_i), 0, cos(theta_j));
2498
2499         G4CutTubs* cutTube = new G4CutTubs(
2500             "CutTube_" + std::to_string(i), rMin, rMax, dz, sPhi, dPhi, pLowNorm, pHighNorm);
2501
2502         //G4Tubs* hole = new G4Tubs("Hole_" + std::to_string(i), 0, holeRadius, TubeR + 1 * cm, 0, 360 * deg);
2503
2504         G4double theta = i * 45 * deg;
2505         G4RotationMatrix* rotation = new G4RotationMatrix();
2506         rotation->rotateY(45 * i * deg);
2507
2508         G4RotationMatrix* rotation_hole = new G4RotationMatrix();
2509         rotation_hole->rotateX(90 * deg);
2510
2511         G4ThreeVector position(-distance * cos(theta), yoffset, distance * sin(theta));
2512         G4Transform3D tr = G4Transform3D(*rotation, position);
2513
2514         //G4SubtractionSolid* cutTubeWithHole = new G4SubtractionSolid(
2515         //    "CutTubeWithHole_" + std::to_string(i), cutTube, hole, rotation_hole, G4ThreeVector(0, 0, 0));
2516
2517
2518         multiUnion->AddNode(*cutTube, tr);
2519     }
2520 }
2521
2522
2523 /**
2524 U4SolidMaker::Snake
2525 --------------------
2526
2527 WaterdistributorpipeMaker::MakeWaterDistributorPartIIIUnion
2528
2529 **/
2530
2531
2532 const G4VSolid* U4SolidMaker::WaterDistributorPartIIIUnion(const char* name_) // static
2533 {
2534     //double m_lowerchimney_rotangle = 0. ;
2535
2536     // Create individual solid components first
2537     double flange_r_inner = 66.5*mm;
2538     double flange_r_outer = 125.*mm;
2539     int numZplanes_flange = 4;
2540     G4double zPlane_flange[] = {-50 *mm, -25 *mm , -25 *mm, 0 *mm};
2541     G4double rInner_flange[] = {0 , 0, 0, 0};
2542     G4double rOuter_flange[] = {flange_r_inner, flange_r_inner, flange_r_outer , flange_r_outer};
2543     auto solidFlange = new G4Polycone("sBotChimneySSEnclosure_Flange",
2544                                     0, 360*deg,
2545                                     numZplanes_flange, zPlane_flange, rInner_flange, rOuter_flange);
2546
2547     double dz1 = 529/2 * mm;
2548     double theta1 = 27*deg;
2549     G4ThreeVector pLowNorm1( sin(theta1), 0, -cos(theta1) );
2550     G4ThreeVector pHighNorm1( -sin(theta1), 0, cos(theta1) );
2551     auto cutTube1 = new G4CutTubs("CutTube1", 0, 66.5 * cos(theta1), dz1, 0, 360*deg, pLowNorm1, pHighNorm1);
2552
2553     double dz2 = 1125.25/2 * mm;
2554     double theta2 = -16*deg;
2555     G4ThreeVector pLowNorm2( sin(theta2), 0, -cos(theta2) );
2556     G4ThreeVector pHighNorm2( -sin(theta2), 0, cos(theta2) );
2557     auto cutTube2 = new G4CutTubs("CutTube2", 0, 66.5 * cos(theta2), dz2, 0, 360*deg, pLowNorm2, pHighNorm2);
2558
2559     double r_tube = 66.5 * mm;
2560     double r_main1 = 283 * mm;
2561     double phi_start1 = 180 * deg;
2562     double phi_delta1 = 135 * deg;
2563     auto torus1 = new G4Torus("torus1", 0, r_tube, r_main1, phi_start1, phi_delta1);
2564
2565     double r_main2 = 352 * mm;
2566     double phi_start2 = 45 * deg;
2567     double phi_delta2 = 45 * deg;
2568     auto torus2 = new G4Torus("torus2", 0, r_tube, r_main2, phi_start2, phi_delta2);
2569
2570     // Component definition arrays
2571     G4VSolid* solids[] = {solidFlange, cutTube1, cutTube2, torus1, torus2};
2572
2573     // Rotation angles around X and Y axes
2574     double rotX_angles[] = {0*deg, 0*deg, 0*deg, -90*deg, 90*deg};
2575     double rotY_angles[] = {0*deg, 27*deg, -16*deg, 180*deg, 0*deg};
2576
2577     // Pre-calculate offsets for CutTubes
2578     double offset1 = dz1 * cos(27*deg);
2579     double offset2 = dz2 * cos(-16*deg);
2580
2581     // Position vectors
2582     G4ThreeVector positions[] = {
2583         G4ThreeVector(0, 0, -19120),                                                          // Flange
2584         G4ThreeVector(-120, 0, -19120 - 50*mm - offset1),                                    // CutTube1
2585         G4ThreeVector(-240 + dz2 * sin(16*deg), 0, -19120 - 50*mm - offset1 * 2 - offset2), // CutTube2
2586         G4ThreeVector(70 - r_main1, 0, -20724),                                             // Torus1
2587         G4ThreeVector(70 - 449 - 283, 0, -21181)                                            // Torus2
2588     };
2589
2590     // Component names for debugging
2591     //const char* component_names[] = {"Flange", "CutTube1", "CutTube2", "Torus1", "Torus2"};
2592
2593     // Create MultiUnion and add all components in a loop
2594     G4MultiUnion* unionSolid = new G4MultiUnion("WaterDistributorPartIIIUnion");
2595
2596     for (int i = 0; i < 5; ++i) {
2597         // Create rotation matrix for this component
2598         G4RotationMatrix rotation;
2599         if (rotX_angles[i] != 0) {
2600             rotation.rotateX(rotX_angles[i]);
2601         }
2602         if (rotY_angles[i] != 0) {
2603             rotation.rotateY(rotY_angles[i]);
2604         }
2605
2606         // Create transform and add node
2607         G4Transform3D transform(rotation, positions[i]);
2608         unionSolid->AddNode(*solids[i], transform);
2609     }
2610
2611     unionSolid->Voxelize();
2612
2613     // Apply final rotation for the entire assembly
2614     //G4RotationMatrix finalRotation;
2615     //finalRotation.rotateZ(-m_lowerchimney_rotangle);
2616
2617     return unionSolid ;
2618 }
2619
2620
2621
2622 /**
2623 U4SolidMaker::MakeLowerWaterDistributorCurvedCutTubes
2624 -------------------------------------------------------
2625
2626 After WaterdistributorpipeMaker::MakeLowerWaterDistributorCurvedCutTubes
2627
2628
2629 2025-12-17 14:14:52.868 INFO  [2884720] [U4SolidMaker::OuterReflectorOrbSubtraction@2690] m_radOuterReflector 20476
2630 2025-12-17 14:14:52.868 INFO  [2884720] [U4SolidMaker::MakeLowerWaterDistributorCurvedCutTubes@2653]
2631                                                     offset1    235.226
2632                                                     offset2    540.830
2633                                         offset1*2 + offset2   1011.281
2634                                       offset1*2 + offset2*2   1552.111
2635  trans_cutTube1 [  -120.000,     0.000,-19405.226]
2636  trans_cutTube2 [   -84.920,     0.000,-20181.281]
2637
2638
2639
2640 ::
2641
2642
2643     In [4]: z = np.array( [-20476, -20181, -19405, -19170] )
2644
2645     In [5]: z - z[0]
2646     Out[5]: array([   0,  295, 1071, 1306])
2647
2648     In [6]: z - z[-1]
2649     Out[6]: array([-1306, -1011,  -235,     0])
2650
2651
2652     20476 - 1552 = 18924
2653
2654     20476 - 235 - 541 = 19700
2655
2656
2657
2658 Start centered::
2659
2660
2661                 +-----+  -
2662                /     /
2663               /     /    offset1
2664              /     /
2665             /  +  /     -                   z = 0
2666            /     /
2667           /     /       offset1
2668          /     /
2669         +-----+         -
2670
2671
2672 Offset by  -19120 - 50*mm = -19170 puts middle at -19170::
2673
2674
2675                 +-----+  -
2676                /     /
2677               /     /    offset1
2678              /     /
2679             /  +  /     - -  - - - - - -    z = -19170
2680            /     /
2681           /     /       offset1
2682          /     /
2683         +-----+         -
2684
2685
2686 Then offset more, -19170 - offset1 puts top at -19170::
2687
2688
2689
2690                             +-----+  -   - - - - - - - - - -   z = -19170
2691                            /     /
2692                           /     /    offset1
2693                          /     /
2694                         /  +  /     - -  - - - - - -    z = -19170 - offset1              = -19170 - 235.22              = -19405.22
2695                        /     /
2696                       /     /       offset1
2697                      /     /
2698                     +-----+         -                   z =  -19170 - 2*offset1          = -19170 - 2*235.22             = -19640.44
2699                      \     \
2700                       \     \
2701                        \     \
2702                         \     \
2703                          \     \
2704                           \  +  \                       z = -19170 - 2*offset1 - offset2 = -19170 - 2*235.22 - 540.830  =  -20181.27
2705                            \     \
2706                             \     \
2707                              \     \                                                                                -------- -20476 --
2708                               \     \
2709                                \     \
2710                                 \     \
2711                                  +-----+                z = -19170 - 2*offset2 - 2*offset2 = -19170 - 2*235.22 - 2*540.830 = -20722.1
2712
2713
2714 **/
2715
2716
2717
2718 std::vector<PlacedSolid> U4SolidMaker::MakeLowerWaterDistributorCurvedCutTubes(double m_lowerchimney_rotangle) {
2719     std::vector<PlacedSolid> cutTubes;
2720
2721     double UNCOINCIDE_MM = ssys::getenvdouble(U4SolidMaker__MakeLowerWaterDistributorCurvedCutTubes_UNCOINCIDE_MM, 0 );
2722
2723     double rMin = 0;
2724     //double rMax = 50*mm;
2725     double sPhi = 0;
2726     double dPhi = 360*deg;
2727
2728
2729     // First section parameters
2730     double dz1 = 529/2 * mm;
2731     double theta1 = 27*deg;
2732     G4ThreeVector pLowNorm1( sin(theta1), 0, -cos(theta1) );
2733     G4ThreeVector pHighNorm1( -sin(theta1), 0, cos(theta1) );
2734     auto cutTube1 = new G4CutTubs("CutTube1", rMin, 66.5 * cos(theta1), dz1, sPhi, dPhi, pLowNorm1, pHighNorm1);
2735
2736     // Second section parameters
2737     double dz2 = 1125.25/2 * mm;
2738     double theta2 = -16*deg;
2739     G4ThreeVector pLowNorm2( sin(theta2), 0, -cos(theta2) );
2740     G4ThreeVector pHighNorm2( -sin(theta2), 0, cos(theta2) );
2741     auto cutTube2 = new G4CutTubs("CutTube2", rMin, 66.5 * cos(theta2), dz2, sPhi, dPhi, pLowNorm2, pHighNorm2);
2742
2743     double offset1 = dz1 * cos(27*deg);
2744     double offset2 = dz2 * cos(-16*deg);
2745
2746     double cos_rot = cos(m_lowerchimney_rotangle);
2747     double sin_rot = sin(m_lowerchimney_rotangle);
2748     double sin_theta2 = sin(16*deg);
2749
2750     // Rotation and translation for cutTube1
2751     G4RotationMatrix rot_cutTube1;
2752     rot_cutTube1.rotateY(27*deg);
2753     rot_cutTube1.rotateZ(-m_lowerchimney_rotangle);
2754
2755     G4ThreeVector trans_cutTube1(-120 * cos_rot, 120 * sin_rot, -19120 - 50*mm - offset1);
2756
2757     cutTubes.push_back({
2758         cutTube1,
2759         trans_cutTube1,
2760         rot_cutTube1,
2761         0,
2762         "lCutTube1",
2763         "pCutTube1"
2764     });
2765
2766     // Rotation and translation for cutTube2
2767     G4RotationMatrix rot_cutTube2;
2768     rot_cutTube2.rotateY(-16*deg);
2769     rot_cutTube2.rotateZ(-m_lowerchimney_rotangle);
2770
2771     G4ThreeVector trans_cutTube2(
2772         2 * (-120 * cos_rot) + dz2 * sin_theta2 * cos_rot,
2773         2 * (120 * sin_rot) - dz2 * sin_theta2 * sin_rot,
2774         -19120 - 50*mm - offset1 * 2 - offset2 + UNCOINCIDE_MM*mm
2775     );
2776
2777
2778     cutTubes.push_back({
2779         cutTube2,
2780         trans_cutTube2,
2781         rot_cutTube2,
2782         1,
2783         "lCutTube2",
2784         "pCutTube2"
2785     });
2786
2787
2788     std::array<double,3> tr1 = {trans_cutTube1.x(), trans_cutTube1.y(), trans_cutTube1.z() };
2789     std::array<double,3> tr2 = {trans_cutTube2.x(), trans_cutTube2.y(), trans_cutTube2.z() };
2790
2791     LOG(info)
2792        << "\n"
2793        << "[" << U4SolidMaker__MakeLowerWaterDistributorCurvedCutTubes_UNCOINCIDE_MM << "]"
2794        << "\n"
2795        << std::setw(60) << " UNCOINCIDE_MM " << std::fixed << std::setw(10) << std::setprecision(3) << UNCOINCIDE_MM
2796        << "\n"
2797        << std::setw(60) << " offset1 " << std::fixed << std::setw(10) << std::setprecision(3) << offset1
2798        << "\n"
2799        << std::setw(60) << " offset2 "  << std::fixed << std::setw(10) << std::setprecision(3) << offset2
2800        << "\n"
2801        << std::setw(60) <<  " offset1*2 + offset2 " << std::fixed << std::setw(10) << std::setprecision(3) << (offset1*2.+offset2)
2802        << "\n"
2803        << std::setw(60) <<  " offset1*2 + offset2*2 " << std::fixed << std::setw(10) << std::setprecision(3) << (offset1*2.+offset2*2.)
2804        << "\n"
2805        << " trans_cutTube1 " << s_bb::Desc_<double,3>(tr1.data())
2806        << "\n"
2807        << " trans_cutTube2 " << s_bb::Desc_<double,3>(tr2.data())
2808        ;
2809
2810     return cutTubes;
2811 }
2812
2813
2814 /**
2815 U4SolidMaker::OuterReflectorOrbSubtraction
2816 --------------------------------------------
2817
2818 Checking::
2819
2820     DetSim1Construction::helper_subtract_water_distributor
2821
2822 **/
2823
2824 const G4VSolid* U4SolidMaker::OuterReflectorOrbSubtraction(const char* name_) // static
2825 {
2826     G4String name = name_ ;
2827
2828     double m_radOuterWater = 20.50*m - 26*mm; // due to back of veto pmt mask is -25.05 mm (See R12860OnlyFrontMaskManager).
2829     double m_thicknessReflector = 2*mm;
2830     double m_radOuterReflector = m_radOuterWater + m_thicknessReflector;
2831
2832
2833     double radOuterReflector = ssys::getenvfloat(U4SolidMaker__OuterReflectorOrbSubtraction_radOuterReflector, m_radOuterReflector );
2834     G4VSolid* solid_output = new G4Orb(name_,  radOuterReflector );
2835
2836
2837     double m_lowerchimney_rotangle = 0. ;
2838
2839   // ==== New: Subtract two inclined cutTubes of lower water distributor ====
2840     auto cutTubes = MakeLowerWaterDistributorCurvedCutTubes(m_lowerchimney_rotangle);
2841
2842     for (size_t i = 0; i < cutTubes.size(); ++i) {
2843         const auto& tube = cutTubes[i];
2844         G4Transform3D transform(tube.rot, tube.pos);
2845
2846         solid_output = new G4SubtractionSolid(
2847             name + "_cutTube" + std::to_string(i+1),
2848             solid_output,
2849             tube.solid,
2850             transform
2851         );
2852     }
2853
2854     G4String solid_output_name = solid_output->GetName();
2855
2856     LOG(info)
2857          << "\n"
2858          << "m_radOuterReflector " << m_radOuterReflector     // 20476 mm
2859          << "\n"
2860          << " [" << U4SolidMaker__OuterReflectorOrbSubtraction_radOuterReflector << "]"
2861          << "\n"
2862          << " radOuterReflector " << radOuterReflector
2863          << "\n"
2864          << " solid_output_name " << solid_output_name
2865          ;
2866
2867     return solid_output ;
2868 }
2869
2870
2871
2872
2873
2874
2875
2876
2877
2878
2879
2880
2881
2882
2883
2884 struct R12860_PMTSolid_Maker
2885 {
2886     R12860_PMTSolid_Maker(double R1, double R1z, double R2, double R3, double H, double H3);
2887
2888     G4VSolid* GetSolid(G4String solidname, double thickness=0.0, char X='\0' );
2889
2890     double m_R1;
2891     double m_R1z;
2892     double m_R1p;
2893     double m_R2;
2894     double m_R3;
2895     double m_theta;
2896     double m_hypotenuse ;
2897     double m_H;
2898     double m_H_1_2;
2899
2900     double m_H3;
2901
2902     double m_delta_torlerance;
2903
2904 };
2905
2906
2907
2908
2909 R12860_PMTSolid_Maker::R12860_PMTSolid_Maker(
2910     double R1,
2911     double R1z,
2912     double R2,
2913     double R3,
2914     double H,
2915     double H3
2916     )
2917     :
2918     m_R1(R1),    // 101.0
2919     m_R1z(R1z),  //  75.5
2920     m_R2(R2),    //  23.0   radius of torus circle
2921     m_R3(R3),    //  42.5
2922     m_H(H),      // 220.0
2923     m_H3(H3)     //  62.0
2924 {
2925     m_H_1_2 = m_H - m_R1z - m_H3;              //    220 - 75.5 - 62 =  82.5
2926     m_theta = atan((m_R2+m_R3)/(m_H_1_2));     //    23+42.5 = 65.5   65.6/82.5 = 0.7951 atan(0.7951)
2927     m_hypotenuse = sqrt(pow(m_H_1_2,2)+pow(m_R2+m_R3,2)) ;
2928     m_R1p = m_hypotenuse - m_R2;  // hypot distance from ellipse center to torus circle
2929
2930     m_delta_torlerance = 1E-2*mm;
2931
2932
2933 }
2934
2935
2936 /**
2937                   ........+.........                          +                +
2938             ..''''        |         ''''..
2939         ..''              |               ''..
2940       .'                  |                   '.            m_R1z
2941     .'                   r1zt                   '.
2942    '                      |                       '
2943   '                       |                        '
2944  '                        +----------r1t----------- +         +
2945   '                      /:\                       '
2946    '                    / : \                     '
2947     '.                 /  :  \                  .'
2948       '.              /   :   \               .'             m_H_1_2           m_H
2949         ''..         /    :    \          ..''
2950             ''''.   /     :     \    .''''
2951                   ./      :      \  .   .
2952                   /               \
2953                  /                 \
2954                 /                   \
2955                +    +     +     +    +
2956                     |           |
2957                     |           |
2958                     |           |
2959                     |           |
2960                     |           |
2961                     |           |
2962                     |           |
2963                     +-----------+
2964      Z
2965      |
2966      +--X
2967
2968
2969  symbol : pmttube_solid_tube(scarf)
2970  name   : R12860_PMTSolidU_pmt_solid_2_Tube
2971  r4t    : 51.1993
2972  h2t/2  : 9.00616
2973  symbol  : pmttube_solid_torus(scarfSub)
2974  name    : R12860_PMTSolidU_pmt_solid_2_Torus
2975  r2t     : 22.999
2976  r2t+r3t : 65.5
2977  symbol  : pmttube_solid_part2
2978  name    : R12860_PMTSolidU_pmt_solid_part2
2979  -h2t/2  : -9.00616(neck-torus-dz)
2980  symbol  : pmttube_solid_1_2(bulb+neck)
2981  name    : R12860_PMTSolidU_pmt_solid_1_2
2982  neck_dz : -73.4938
2983  total_torus_dz : -82.5
2984  symbol  : pmttube_solid_1_2_3((bulb+neck)+endtube)
2985  name    : R12860_PMTSolidU_pmt_solid
2986  endtube_dz : -113.501
2987
2988
2989 POINT=-65.5,0,-82.5,65.5,0,-82.5 BBOX=-51.19,0,-82.4999,51.19,0,-64.48764,-51.19,0,-64.48764,51.19,0,-82.4999   KLUDGE=1 cxt_min.sh pdb
2990
2991
2992 **/
2993
2994
2995
2996
2997 G4VSolid* R12860_PMTSolid_Maker::GetSolid(G4String solidname, double thickness, char X)
2998 {
2999     // Calculate Parameter first
3000     double r1t = m_R1 + thickness;
3001     double r1zt = m_R1z + thickness;
3002     double r2t = m_R2 - thickness;
3003     double r3t = m_R3 + thickness;
3004     double r1zp = (m_R1p+thickness);
3005     double r4t = r1zp * sin(m_theta);
3006
3007     double h1t = r1zp * cos(m_theta);   // above necktop height (up to ellipse equator)
3008     double h2t = r2t * cos(m_theta);    // below necktop height (down to torus plane line)
3009
3010     double h3t = m_H3 + thickness;
3011
3012
3013     double scarf_radius = r4t+m_delta_torlerance ;
3014     double scarf_halfheight = h2t/2+m_delta_torlerance ;
3015
3016     double endtube_radius = r3t+m_delta_torlerance ;
3017     double endtube_halfheight = h3t/2+m_delta_torlerance ;
3018
3019
3020     std::cout
3021         << " thickness           : " << std::setw(10) << thickness          << "\n"
3022         << " m_delta_torlerance    : " << std::setw(10) << m_delta_torlerance   << "\n"
3023         << " m_R1                : " << std::setw(10) << m_R1               << " (input ellipse width)"                  << "\n"
3024         << " r1t                 : " << std::setw(10) << r1t                << " (m_R1 + thickness)"                     << "\n"
3025         << " m_R1z               : " << std::setw(10) << m_R1z              << " (input ellipse height)"                 << "\n"
3026         << " r1zt                : " << std::setw(10) << r1zt               << " (m_R1z + thickness)"                    << "\n"
3027         << " m_R2                : " << std::setw(10) << m_R2               << " (input torus circle radius)"            << "\n"
3028         << " m_R3                : " << std::setw(10) << m_R3               << " (input endtube radius)"                 << "\n"
3029         << " m_H                 : " << std::setw(10) << m_H                << " (input total height)"                   << "\n"
3030         << " m_H3                : " << std::setw(10) << m_H3               << " (input endtube height)"                 << "\n"
3031         << " m_H_1_2             : " << std::setw(10) << m_H_1_2            << " (m_H - m_R1z - m_H3)"                   << "\n"
3032         << " m_theta             : " << std::setw(10) << m_theta            << " (ellipse ^ torus centers)"              << "\n"
3033         << " m_hypotenuse        : " << std::setw(10) << m_hypotenuse       << " (dist from ellipse to torus center)"    << "\n"
3034         << " m_R1p               : " << std::setw(10) << m_R1p              << " (m_hypotenuse - m_R2)"                  << "\n"
3035         << " r1zp                : " << std::setw(10) << r1zp               << " (m_R1p + thickness)"                    << "\n"
3036         << " r4t                 : " << std::setw(10) << r4t                << " (r1zp*sin(m_theta))"                    << "\n"
3037         << " scarf_radius        : " << std::setw(10) << scarf_radius       << " (r4t+m_delta_torlerance)"                 << "\n"
3038         << " r2t                 : " << std::setw(10) << r2t                << " (m_R2 - thickness)"                     << "\n"
3039         << " h1t                 : " << std::setw(10) << h1t                << " (r1zp*cos(m_theta)) above neck height " << "\n"
3040         << " h2t                 : " << std::setw(10) << h2t                << " (r2t*cos(m_theta))  neck h"             << "\n"
3041         << " scarf_halfheight    : " << std::setw(10) << scarf_halfheight   << " (h2t/2+m_delta_torlerance)"               << "\n"
3042         << " endtube_radius      : " << std::setw(10) << endtube_radius     << " (r3t+m_delta_torlerance)"                 << "\n"
3043         << " h3t                 : " << std::setw(10) << h3t                << " (m_H3 + thickness)"                     << "\n"
3044         << " endtube_halfheight  : " << std::setw(10) << endtube_halfheight << " (h3t/2+m_delta_torlerance)"               << "\n"
3045         ;
3046
3047
3048
3049     // Show variables
3050     G4cout << "r1t: " << r1t/mm << " mm" << G4endl;
3051     G4cout << "r4t: " << r4t/mm << " mm" << G4endl;
3052     G4cout << "r2t: " << r2t/mm << " mm" << G4endl;
3053     G4cout << "r3t: " << r3t/mm << " mm" << G4endl;
3054     G4cout << "h1t: " << h1t/mm << " mm" << G4endl;
3055     G4cout << "h2t: " << h2t/mm << " mm" << G4endl;
3056     G4cout << "h3t: " << h3t/mm << " mm" << G4endl;
3057
3058     // Construct the PMT Solid
3059     // * PART 1
3060     // G4Sphere* pmttube_solid_sphere = new G4Sphere(
3061     //                                         solidname+"_1_Sphere",
3062     //                                         0*mm, // R min
3063     //                                         r1t, // R max
3064     //                                         0*deg, // Start Phi
3065     //                                         360*deg, // Delta Phi
3066     //                                         0*deg, // Start Theta
3067     //                                         180*deg  // Delta Theta
3068     //                                         );
3069     G4Ellipsoid* pmttube_solid_sphere = new G4Ellipsoid(
3070                                             solidname+"_1_Ellipsoid",
3071                                             r1t, // pxSemiAxis
3072                                             r1t, // pySemiAxis
3073                                             r1zt // pzSemiAxis
3074                                             );
3075     // * PART 2  : neck
3076
3077
3078     G4Tubs* pmttube_solid_tube = new G4Tubs(
3079                                     solidname+"_2_Tube",
3080                                     0*mm,  /* inner */
3081                                     scarf_radius, /* pmt_r */
3082                                     scarf_halfheight, /* part 2 h */
3083                                     0*deg,
3084                                     360*deg);
3085
3086     std::cout
3087          << " symbol : pmttube_solid_tube(scarf)" << "\n"
3088          << " name   : " << pmttube_solid_tube->GetName() << "\n"
3089          << " r4t    : " << r4t << "\n"
3090          << " h2t/2  : " << h2t/2 << "\n"
3091          ;
3092
3093     G4Torus* pmttube_solid_torus = new G4Torus(
3094                                         solidname+"_2_Torus",
3095                                         0*mm,  // R min
3096                                         r2t+m_delta_torlerance, // R max
3097                                         (r2t+r3t), // Swept Radius
3098                                         -0.01*deg,
3099                                         360.01*deg);
3100
3101     std::cout
3102          << " symbol  : pmttube_solid_torus(scarfSub) " << "\n"
3103          << " name    : " << pmttube_solid_torus->GetName() << "\n"
3104          << " r2t     : " << r2t << "\n"
3105          << " r2t+r3t : " << r2t+r3t << "\n"
3106          ;
3107
3108
3109
3110
3111
3112     [[maybe_unused]] G4VSolid* pmttube_solid_part2 = nullptr ;
3113
3114     G4VSolid* neck = nullptr ;
3115
3116     if( X == 'U' )
3117     {
3118         pmttube_solid_part2 = new G4UnionSolid(
3119                                             solidname+"_part2",
3120                                             pmttube_solid_tube,
3121                                             pmttube_solid_torus,
3122                                             0,
3123                                             G4ThreeVector(0,0,-h2t/2)
3124                                             );
3125
3126         neck = pmttube_solid_part2 ;   // debug union : allows to see the position of the torus relative to the rest
3127     }
3128     else if( X == 'K' )
3129     {
3130         neck = pmttube_solid_tube ;  // KLUDGE : DONT SUBTRACT TORUS
3131     }
3132     else if( X == 'P' )
3133     {
3134         G4double phiStart = 0.00*deg ;
3135         G4double phiTotal = 360.00*deg ;
3136         G4int numZPlanes = 2 ;
3137         G4double zPlane[] = { -scarf_halfheight, scarf_halfheight  } ;
3138         G4double rInner[] = {  0.0             , 0.0   } ;
3139         G4double rOuter[] = {  endtube_radius  , scarf_radius } ;
3140
3141         G4Polycone* polycone_neck = new G4Polycone(
3142                                    solidname+"_part2",
3143                                    phiStart,
3144                                    phiTotal,
3145                                    numZPlanes,
3146                                    zPlane,
3147                                    rInner,
3148                                    rOuter
3149                                    );
3150
3151         neck = polycone_neck ;
3152     }
3153     else
3154     {
3155         pmttube_solid_part2 = new G4SubtractionSolid(
3156                                             solidname+"_part2",
3157                                             pmttube_solid_tube,
3158                                             pmttube_solid_torus,
3159                                             0,
3160                                             G4ThreeVector(0,0,-h2t/2)
3161                                             );
3162         neck = pmttube_solid_part2 ;   // original monstrosity
3163     }
3164
3165     std::cout
3166          << " symbol  : pmttube_solid_part2 " << "\n"
3167          << " name    : " << ( pmttube_solid_part2 ? pmttube_solid_part2->GetName() : "-" ) << "\n"
3168          << " -h2t/2  : " << -h2t/2 << "(neck-torus-dz)\n"
3169          ;
3170
3171
3172     // * PART 3  : EndTube
3173
3174
3175
3176     G4Tubs* pmttube_solid_end_tube = new G4Tubs(
3177                                     solidname+"_3_EndTube",
3178                                     0*mm,  /* inner */
3179                                     endtube_radius, //21*cm/2, /* pmt_r */
3180                                     endtube_halfheight, //30*cm/2, /* pmt_h */
3181                                     0*deg,
3182                                     360*deg);
3183
3184
3185
3186     // * PART 1 + 2   : bulb + neck
3187     G4UnionSolid* pmttube_solid_1_2 = new G4UnionSolid(
3188                                             solidname+"_1_2",
3189                                             pmttube_solid_sphere,
3190                                             neck,
3191                                             0,
3192                                             G4ThreeVector(0, 0, -(h1t+h2t/2))
3193                                             );
3194
3195      std::cout
3196           << " symbol  : pmttube_solid_1_2(bulb+neck) \n"
3197           << " name    : " << pmttube_solid_1_2->GetName() << "\n"
3198           << " neck_dz : " << -(h1t+h2t/2) << "\n"
3199           << " total_torus_dz : " << -h2t/2 + -(h1t+h2t/2) << "\n"
3200           ;
3201
3202     // * PART 1+2 + 3
3203     [[maybe_unused]] G4UnionSolid* pmttube_solid_1_2_3 = new G4UnionSolid(
3204                                             solidname,
3205                                             pmttube_solid_1_2,
3206                                             pmttube_solid_end_tube,
3207                                             0,
3208                                             G4ThreeVector(0,0,
3209                                                 -(m_H_1_2+h3t*0.50))
3210                                             );
3211
3212      std::cout
3213           << " symbol  : pmttube_solid_1_2_3((bulb+neck)+endtube) \n"
3214           << " name    : " << pmttube_solid_1_2_3->GetName() << "\n"
3215           << " endtube_dz : " << -(m_H_1_2+h3t*0.50) << "\n"
3216           ;
3217
3218
3219     return pmttube_solid_1_2_3;
3220 }
3221
3222
3223
3224 const G4VSolid* U4SolidMaker::R12860_PMTSolid(const char* name_) // static
3225 {
3226     const char* PREFIX = "R12860_PMTSolid" ;
3227     assert( sstr::StartsWith(name_,PREFIX ));
3228     const char* suffix = name_ + strlen(PREFIX);
3229     G4String solidname = name_ ;
3230
3231     // Ham8inchPMTManager::init_variables
3232
3233     double m_pmt_r, m_pmt_h, m_z_equator ;
3234
3235     m_pmt_r = 101.*mm;
3236     m_pmt_h = 220.*mm;
3237     m_z_equator = 75.5*mm; // From top to equator
3238
3239     R12860_PMTSolid_Maker* m_pmtsolid_maker = new R12860_PMTSolid_Maker(
3240                             m_pmt_r,       // m_R1     101.0
3241                             m_z_equator,   // m_R1z     75.5
3242                             23*mm,         // m_R2      23.0
3243                             42.5*mm,       // m_R3      42.5
3244                             m_pmt_h,       // m_H       220.0
3245                             62.*mm         // m_H3      62.0
3246                             );
3247
3248
3249
3250
3251     double thickness = 1E-3*mm ;
3252
3253     bool with_suffix = suffix && ( suffix[0] == 'U' || suffix[0] == 'K' || suffix[0] == 'P' ) ;
3254
3255     return m_pmtsolid_maker->GetSolid(solidname + "_pmt_solid", thickness, ( with_suffix ? suffix[0] : '\0' ) );
3256 }
3257
3258
3259
3260
3261
3262
3263
3264
3265
3266
3267 struct EMFCoil
3268 {
3269     static const double EMF_RC_32[32] ;
3270     static G4Polycone* BuildLProfileSectorPolycone( const std::string& solidName, double Rout_mm, double L_mm, double t_mm, double phiStart, double phiTotal, bool flangeUp = false );
3271     static const G4VSolid* Example(int i);
3272 };
3273
3274 const double EMFCoil::EMF_RC_32[32] = {
3275         3.843,   6.509,   9.185,  11.029,
3276        13.210,  14.924,  16.325,  17.495,
3277        18.465,  19.284,  19.953,  20.502,
3278        20.931,  21.241,  21.450,  21.559,
3279        21.559,  21.450,  21.240,  20.930,
3280        20.501,  19.951,  19.282,  18.463,
3281        17.495,  16.325,  14.924,  13.210,
3282        11.053,   9.185,   6.509,   3.843
3283     };
3284
3285
3286
3287 G4Polycone* EMFCoil::BuildLProfileSectorPolycone(
3288     const std::string& solidName,
3289     double Rout_mm,
3290     double L_mm,
3291     double t_mm,
3292     double phiStart,
3293     double phiTotal,
3294     bool flangeUp
3295 ){
3296
3297     const double H     = L_mm;
3298
3299     const int    NZ    = 4;
3300
3301     double zPlane[NZ];
3302     double rMin[NZ];
3303     double rMax[NZ];
3304
3305     if (!flangeUp) {
3306
3307         zPlane[0] = -H*0.5;
3308         zPlane[1] = -H*0.5 + t_mm;
3309         zPlane[2] = -H*0.5 + t_mm;
3310         zPlane[3] = +H*0.5;
3311
3312         // flange：r ∈ [Rout - L, Rout]
3313         rMin[0] = Rout_mm - L_mm;
3314         rMin[1] = Rout_mm - L_mm;
3315
3316         // web：r ∈ [Rout - t, Rout]
3317         rMin[2] = Rout_mm - t_mm;
3318         rMin[3] = Rout_mm - t_mm;
3319     } else {
3320
3321         zPlane[0] = -H*0.5;
3322         zPlane[1] = -H*0.5 + (H - t_mm);
3323         zPlane[2] = -H*0.5 + (H - t_mm);
3324         zPlane[3] = +H*0.5;
3325
3326         rMin[0] = Rout_mm - t_mm;
3327         rMin[1] = Rout_mm - t_mm;
3328         rMin[2] = Rout_mm - L_mm;
3329         rMin[3] = Rout_mm - L_mm;
3330     }
3331
3332
3333     for (int i = 0; i < NZ; ++i) rMax[i] = Rout_mm;
3334
3335
3336     for (int i = 0; i < NZ; ++i)
3337     {
3338         std::cout
3339            << "EMFCoil::BuildLProfileSectorPolycone"
3340            << " i " << std::setw(2) << i
3341            << " z    " << std::setw(10) << zPlane[i]
3342            << " rMin " << std::setw(10) << rMin[i]
3343            << " rMax " << std::setw(10) << rMax[i]
3344            << "\n"
3345            ;
3346     }
3347
3348     return new G4Polycone(
3349         solidName.c_str(),
3350         phiStart * rad,
3351         phiTotal * rad,
3352         NZ, zPlane, rMin, rMax
3353     );
3354 }
3355
3356
3357
3358 const G4VSolid* EMFCoil::Example(int i) // static
3359 {
3360     // WaterPoolConstruction::makeEMFCoilHolders_Polycone
3361     const double L_leg_mm  = 56.0;
3362     const double t_mm      = 4.0;
3363     const double eps_r_mm  = 2.0;
3364
3365     const double R0_m = EMF_RC_32[i] + eps_r_mm * 1e-3;
3366
3367     double phiS = 89.60*deg ;
3368     double phiE = 134.72*deg ;
3369
3370     const double PHI_EPS = 0.01 * deg;
3371     const double s = phiS + PHI_EPS;
3372     const double e = phiE - PHI_EPS;
3373     const double dphi = e - s;
3374
3375     std::string solidName = "example" ;
3376     double Rout_mm = R0_m * 1000.0 ;
3377     double L_mm = L_leg_mm ;
3378     double phiStart = s ;
3379     double phiTotal = dphi ;
3380     bool flangeUp = false ;
3381
3382     return BuildLProfileSectorPolycone(solidName, Rout_mm, L_mm, t_mm, phiStart, phiTotal, flangeUp );
3383 }
3384
3385
3386 const G4VSolid* U4SolidMaker::LProfileSectorPolycone( const char* )
3387 {
3388     return EMFCoil::Example(16);
3389 }