diff --git a/GeometryService/inc/ProductionTargetMaker.hh b/GeometryService/inc/ProductionTargetMaker.hh index 97eeb9d73e..2f1bd98cf6 100644 --- a/GeometryService/inc/ProductionTargetMaker.hh +++ b/GeometryService/inc/ProductionTargetMaker.hh @@ -6,6 +6,8 @@ #include #include #include +#include +#include namespace mu2e { class SimpleConfig; } namespace mu2e { class ProductionTarget; } @@ -14,15 +16,20 @@ namespace mu2e { class ProductionTargetMaker { public: + // Registry map for model makers — allows scaling without modifying make() + using MakerFunction = std::function(const SimpleConfig&, double)>; + static const std::map& getMakerRegistry(); static std::unique_ptr make(const SimpleConfig& config, double solenoidOffset); static const int tier1{1}; // version 2 is the low density hayman static const int hayman_v_2_0{3}; + static const int stickman_v_1_0{4}; static std::unique_ptr makeTier1(const SimpleConfig& config, double solenoidOffset); static std::unique_ptr makeHayman_v_2_0(const SimpleConfig& config, double solenoidOffset); + static std::unique_ptr makeStickman_v_1_0(const SimpleConfig& config, double solenoidOffset); diff --git a/GeometryService/src/GeometryService.cc b/GeometryService/src/GeometryService.cc index 5edbe1145a..d8dd55a1c9 100644 --- a/GeometryService/src/GeometryService.cc +++ b/GeometryService/src/GeometryService.cc @@ -238,17 +238,8 @@ namespace mu2e { addDetector(PSVacuumMaker::make(*_config, ps, pse, vacPS_TS_z)); - //addDetector(PSShieldMaker::make(*_config, ps.psEndRefPoint(), prodTarget.position())); - - if (_config->getString("targetPS_model") == "MDC2018"){ - // std::cout << "adding Tier1 in GeometryService" << std::endl; - addDetector(PSShieldMaker::make(*_config, ps.psEndRefPoint(), prodTarget.position())); - } else - if (_config->getString("targetPS_model") == "Hayman_v_2_0"){ - // std::cout << " adding Hayman in GeometryService" << std::endl; - addDetector(PSShieldMaker::make(*_config, ps.psEndRefPoint(), prodTarget.haymanProdTargetPosition())); - } else - {throw cet::exception("GEOM") << " " << static_cast(__func__) << " illegal production target version specified in GeometryService_service = " << _config->getString("targetPS_model") << std::endl;} + // Use ProductionTarget's built-in method to get the correct position based on model type + addDetector(PSShieldMaker::make(*_config, ps.psEndRefPoint(), prodTarget.targetPositionByVersion())); diff --git a/GeometryService/src/ProductionTargetMaker.cc b/GeometryService/src/ProductionTargetMaker.cc index 81d1f5f701..fbc1161300 100644 --- a/GeometryService/src/ProductionTargetMaker.cc +++ b/GeometryService/src/ProductionTargetMaker.cc @@ -8,20 +8,27 @@ namespace mu2e { - std::unique_ptr ProductionTargetMaker::make(const SimpleConfig& c, double solenoidOffset) { - - - if (c.getString("targetPS_model") == "MDC2018"){ - // std::cout << "making Tier1 in maker" << std::endl; - return makeTier1(c, solenoidOffset); - } else - if (c.getString("targetPS_model") == "Hayman_v_2_0"){ - // std::cout << " making Hayman in Maker" << std::endl; - return makeHayman_v_2_0(c, solenoidOffset); - } else - {throw cet::exception("GEOM") << " illegal production target version specified = " << c.getInt("targetPS_version") << std::endl;} - return 0; + // Registry of target model makers + const std::map& ProductionTargetMaker::getMakerRegistry() { + static const std::map registry = { + {"MDC2018", [](const SimpleConfig& c, double offset) { return makeTier1(c, offset); }}, + {"Hayman_v_2_0", [](const SimpleConfig& c, double offset) { return makeHayman_v_2_0(c, offset); }}, + {"Stickman_v_1_0", [](const SimpleConfig& c, double offset) { return makeStickman_v_1_0(c, offset); }} + }; + return registry; + } + std::unique_ptr ProductionTargetMaker::make(const SimpleConfig& c, double solenoidOffset) { + const std::string modelName = c.getString("targetPS_model"); + const auto& registry = getMakerRegistry(); + + auto it = registry.find(modelName); + if (it != registry.end()) { + return it->second(c, solenoidOffset); + } + + throw cet::exception("GEOM") + << "illegal production target model specified = " << modelName << std::endl; } std::unique_ptr ProductionTargetMaker::makeTier1(const SimpleConfig& c, double solenoidOffset){ @@ -390,5 +397,197 @@ namespace mu2e { return tgtPS; } + std::unique_ptr ProductionTargetMaker::makeStickman_v_1_0(const SimpleConfig& c, double solenoidOffset){ + + // Read the plate and fin parameters + std::vector plateMaterial; + std::vector plateROut; + std::vector plateFinAngles; + std::vector plateThickness; + std::vector plateLugThickness; + + c.getVectorString("targetPS_plateMaterial", plateMaterial); + c.getVectorDouble("targetPS_rOut", plateROut); + c.getVectorDouble("targetPS_plateFinAngles", plateFinAngles); + for_each(plateFinAngles.begin(), plateFinAngles.end(), [](double& elem){elem *= CLHEP::degree;}); + c.getVectorDouble("targetPS_plateThickness", plateThickness); + c.getVectorDouble("targetPS_plateLugThickness", plateLugThickness); + + const int nPlates = c.getInt("targetPS_numberOfPlates"); + const int nFins = c.getInt("targetPS_nStickmanFins"); + + if (plateMaterial.size() != static_cast(nPlates)) { + throw cet::exception("GEOM") + << "targetPS_plateMaterial size mismatch: expected " << nPlates + << ", got " << plateMaterial.size(); + } + if (plateROut.size() != static_cast(nPlates)) { + throw cet::exception("GEOM") + << "targetPS_rOut size mismatch: expected " << nPlates + << ", got " << plateROut.size(); + } + if (plateThickness.size() != static_cast(nPlates)) { + throw cet::exception("GEOM") + << "targetPS_plateThickness size mismatch: expected " << nPlates + << ", got " << plateThickness.size(); + } + if (plateLugThickness.size() != static_cast(nPlates)) { + throw cet::exception("GEOM") + << "targetPS_plateLugThickness size mismatch: expected " << nPlates + << ", got " << plateLugThickness.size(); + } + if (plateFinAngles.size() != static_cast(nFins)) { + throw cet::exception("GEOM") + << "targetPS_plateFinAngles size mismatch: expected " << nFins + << ", got " << plateFinAngles.size(); + } + + // Build parameter structs for Stickman constructor + StickmanEnvelopeParams envelopeParams{ + c.getDouble("targetPS_productionTargetMotherOuterRadius"), + c.getDouble("targetPS_productionTargetMotherHalfLength"), + c.getDouble("targetPS_rotX") * CLHEP::degree, + c.getDouble("targetPS_rotY") * CLHEP::degree, + c.getDouble("targetPS_rotZ") * CLHEP::degree, + c.getDouble("targetPS_halfStickmanLength"), + CLHEP::Hep3Vector(solenoidOffset, + 0, + c.getDouble("productionTarget.zNominal")) + + c.getHep3Vector("productionTarget.offset"), + c.getString("targetPS_targetVacuumMaterial") + }; + + StickmanPlateParams plateParams{ + nPlates, + plateMaterial, + plateROut, + nFins, + plateFinAngles, + c.getDouble("targetPS_plateFinOuterRadius"), + c.getDouble("targetPS_plateFinWidth"), + c.getDouble("targetPS_plateCenterToLugCenter"), + c.getDouble("targetPS_plateLugInnerRadius"), + c.getDouble("targetPS_plateLugOuterRadius"), + plateThickness, + plateLugThickness + }; + + StickmanRodParams rodParams{ + c.getString("targetPS_rodMaterial"), + c.getDouble("targetPS_rodRadius") + }; + + StickmanSpacerParams spacerParams{ + c.getString("targetPS_spacerMaterial"), + c.getDouble("targetPS_spacerHalfLength"), + c.getDouble("targetPS_spacerOuterRadius"), + c.getDouble("targetPS_spacerInnerRadius") + }; + + StickmanSupportRingParams supportRingParams{ + c.getString("targetPS_supportRingMaterial"), + c.getDouble("targetPS_supportRingLength"), + c.getDouble("targetPS_supportRingInnerRadius"), + c.getDouble("targetPS_supportRingOuterRadius"), + c.getDouble("targetPS_supportRingLugOuterRadius"), + c.getDouble("targetPS_supportRingCutoutOffset") + }; + + std::unique_ptr tgtPS + (new ProductionTarget( + c.getString("targetPS_model","NULL"), + c.getInt("targetPS_version"), + envelopeParams, + plateParams, + rodParams, + spacerParams, + supportRingParams + )); + + // Create configuration parameters struct + StickmanConfigParams configParams; + + // Configure plate fillet parameters (only if fillets will be used) + configParams.addFilletToPlateCore = c.getBool("targetPS_addFilletToPlateCore"); + configParams.addFilletToPlateLug = c.getBool("targetPS_addFilletToPlateLug"); + if(configParams.addFilletToPlateCore || configParams.addFilletToPlateLug) { + configParams.plateFilletRadius = c.getDouble("targetPS_plateFilletRadius"); + } + + // Configure support ring lug fillet parameters (only if fillets will be used) + configParams.addFilletToSupportRingLug = c.getBool("targetPS_addFilletToSupportRingLug"); + if(configParams.addFilletToSupportRingLug) { + configParams.supportRingLugFilletRadius = c.getDouble("targetPS_supportRingLugFilletRadius"); + } + + // Configure support ring cutout parameters (only if cutouts will be used) + configParams.addCutoutToSupportRing = c.getBool("targetPS_addCutoutToSupportRing"); + if(configParams.addCutoutToSupportRing) { + configParams.nSupportRingCutouts = c.getInt("targetPS_nSupportRingCutouts"); + c.getVectorDouble("targetPS_supportRingCutoutAngles", configParams.supportRingCutoutAngles); + if (configParams.supportRingCutoutAngles.size() != static_cast(configParams.nSupportRingCutouts)) { + throw cet::exception("GEOM") + << "targetPS_supportRingCutoutAngles size mismatch: expected " + << configParams.nSupportRingCutouts << ", got " << configParams.supportRingCutoutAngles.size(); + } + std::for_each(configParams.supportRingCutoutAngles.begin(), configParams.supportRingCutoutAngles.end(), [](double& elem){elem *= CLHEP::degree;}); + configParams.supportRingCutoutInnerRadius = c.getDouble("targetPS_supportRingCutoutInnerRadius"); + configParams.supportRingCutoutTilt = c.getDouble("targetPS_supportRingCutoutTilt") * CLHEP::degree; + } + + // Configure support wheel/bicycle wheel parameters (reused from Hayman) + configParams.supportsBuild = c.getBool("targetPS.supports.build", false); + if(configParams.supportsBuild) { + //support wheel parameters + configParams.supportWheelRIn = c.getDouble("targetPS.supports.wheel.rIn"); + configParams.supportWheelROut = c.getDouble("targetPS.supports.wheel.rOut"); + configParams.supportWheelHL = c.getDouble("targetPS.supports.wheel.halfLength"); + configParams.supportWheelMaterial = c.getString("targetPS.supports.wheel.material"); + //number of support rods and wires + configParams.nSpokesPerSide = c.getInt("targetPS.supports.nSpokes"); + //features on the wheel near each support rod + c.getVectorDouble("targetPS.supports.features.angles", configParams.supportWheelFeatureAngles, configParams.nSpokesPerSide); + c.getVectorDouble("targetPS.supports.features.arcs" , configParams.supportWheelFeatureArcs , configParams.nSpokesPerSide); + c.getVectorDouble("targetPS.supports.features.rIns" , configParams.supportWheelFeatureRIns , configParams.nSpokesPerSide); + //support wheel rods parameters + c.getVectorDouble("targetPS.supports.rods.halfLength", configParams.supportWheelRodHL, configParams.nSpokesPerSide); + c.getVectorDouble("targetPS.supports.rods.offset", configParams.supportWheelRodOffset, configParams.nSpokesPerSide); + c.getVectorDouble("targetPS.supports.rods.pinOffset", configParams.supportWheelRodPinOffset, configParams.nSpokesPerSide); + c.getVectorDouble("targetPS.supports.rods.radius", configParams.supportWheelRodRadius, configParams.nSpokesPerSide); + c.getVectorDouble("targetPS.supports.rods.radialOffset", configParams.supportWheelRodRadialOffset, configParams.nSpokesPerSide); + c.getVectorDouble("targetPS.supports.rods.wireOffset.downstream", configParams.supportWheelRodWireOffsetD, configParams.nSpokesPerSide); + c.getVectorDouble("targetPS.supports.rods.wireOffset.upstream", configParams.supportWheelRodWireOffsetU, configParams.nSpokesPerSide); + c.getVectorDouble("targetPS.supports.rods.angles", configParams.supportWheelRodAngles, configParams.nSpokesPerSide); + //support wire (spokes) parameters + c.getVectorDouble("targetPS.supports.spokes.targetAngles.downstream", configParams.spokeTargetAnglesD, configParams.nSpokesPerSide); + c.getVectorDouble("targetPS.supports.spokes.targetAngles.upstream", configParams.spokeTargetAnglesU, configParams.nSpokesPerSide); + configParams.spokeRadius = 0.5*c.getDouble("targetPS.supports.spokes.diameter"); + configParams.spokeMaterial = c.getString("targetPS.supports.spokes.material"); + //check that all support wheel vectors are the correct size + if (configParams.supportWheelFeatureAngles.size() != static_cast(configParams.nSpokesPerSide) || + configParams.supportWheelFeatureArcs.size() != static_cast(configParams.nSpokesPerSide) || + configParams.supportWheelFeatureRIns.size() != static_cast(configParams.nSpokesPerSide) || + configParams.supportWheelRodHL.size() != static_cast(configParams.nSpokesPerSide) || + configParams.supportWheelRodOffset.size() != static_cast(configParams.nSpokesPerSide) || + configParams.supportWheelRodPinOffset.size() != static_cast(configParams.nSpokesPerSide) || + configParams.supportWheelRodRadius.size() != static_cast(configParams.nSpokesPerSide) || + configParams.supportWheelRodRadialOffset.size() != static_cast(configParams.nSpokesPerSide) || + configParams.supportWheelRodWireOffsetD.size() != static_cast(configParams.nSpokesPerSide) || + configParams.supportWheelRodWireOffsetU.size() != static_cast(configParams.nSpokesPerSide) || + configParams.supportWheelRodAngles.size() != static_cast(configParams.nSpokesPerSide) || + configParams.spokeTargetAnglesD.size() != static_cast(configParams.nSpokesPerSide) || + configParams.spokeTargetAnglesU.size() != static_cast(configParams.nSpokesPerSide)) { + throw cet::exception("GEOM") + << "Support configuration vector size mismatch for targetPS.supports.nSpokes=" + << configParams.nSpokesPerSide; + } + } + + // Configure the target with the parameters + tgtPS->configureStickman(configParams); + + return tgtPS; + } + } // namespace mu2e diff --git a/Mu2eG4/geom/ProductionTarget_Stickman_v1_0.txt b/Mu2eG4/geom/ProductionTarget_Stickman_v1_0.txt new file mode 100644 index 0000000000..c764852f0a --- /dev/null +++ b/Mu2eG4/geom/ProductionTarget_Stickman_v1_0.txt @@ -0,0 +1,149 @@ +// A geometry file for using the version 1 Stickman Production Target. + +// Support wheel and wire geometry is inherited from the existing Hayman production target geometry, +// others are overhauled to use the new Stickman geometry. + +// The Hayman dependence is intentionally removed. + +int PSStickman.verbosityLevel = 0; //5; + +int targetPS_version = 4; +string targetPS_model = "Stickman_v_1_0"; + +// Nominal target position in the Mu2e coordinate system: at PS center in XY, at the given Z: +double productionTarget.zNominal = -6164.5; +// Optional shift of production target from the nominal +vector productionTarget.offset = { 0., 0., 0.}; + +// Mother volume +double targetPS_productionTargetMotherOuterRadius = 200.; //mm +double targetPS_productionTargetMotherHalfLength = 120.; //mm full target length 232.2 mm, + // allow the mother volume to be slightly larger + +// Overall rotation +double targetPS_rotX = 0.0; +double targetPS_rotY = 14.0; +double targetPS_rotZ = 0.0; +// Dimension +double targetPS_halfStickmanLength = 116.1; // mm. 35 * 6 mm / 2 + 8.1 mm + 3 mm. Actual half length. +// Configurations +bool targetPS.visible = true; +bool targetPS.solid = false; +string targetPS_targetVacuumMaterial = "PSVacuum"; + +// Start with target plates. Chose to specify certain parameters for each plate, so can easily play with plate-by-plate dimensions later on +int targetPS_numberOfPlates = 35; +vector targetPS_plateMaterial = { // counting along mu2e +z direction, i.e. from the proton downstream side to the upstream side + "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 5 + "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 10 + "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 15 + "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 20 + "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 25 + "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718", // 30 + "Inconel718", "Inconel718", "Inconel718", "Inconel718", "Inconel718" // 35 +}; // does not separate fin & core material here +vector targetPS_rOut = { + 3.15, 3.15, 3.15, 3.15, 3.15, // 5 + 3.15, 3.15, 3.15, 3.15, 3.15, // 10 + 3.15, 3.15, 3.15, 3.15, 3.15, // 15 + 3.15, 3.15, 3.15, 3.15, 3.15, // 20 + 3.15, 3.15, 3.15, 3.15, 3.15, // 25 + 3.15, 3.15, 3.15, 3.15, 3.15, // 30 + 3.15, 3.15, 3.15, 3.15, 3.15 // 35 +}; // in mm, core part radius, one value for each +int targetPS_nStickmanFins = 3; // since the plates are installed on rods, the number of fins is always the same +vector targetPS_plateFinAngles = {285.,165.,45.}; // degrees +double targetPS_plateFinOuterRadius = 18.075; // mm, (39.20-3.05)/2, extends to lug-to-center radius - lug inner radius +double targetPS_plateFinWidth = 2.0; // mm, all fin has same width +double targetPS_plateCenterToLugCenter = 19.6; // mm, center of plate to center of lug +double targetPS_plateLugInnerRadius = 1.525; // mm, inner radius of the lug +double targetPS_plateLugOuterRadius = 3.0; // mm, outer radius of the lug +vector targetPS_plateThickness = { + 5.0, 5.0, 5.0, 5.0, 5.0, // 5 + 5.0, 5.0, 5.0, 5.0, 5.0, // 10 + 5.0, 5.0, 5.0, 5.0, 5.0, // 15 + 5.0, 5.0, 5.0, 5.0, 5.0, // 20 + 5.0, 5.0, 5.0, 5.0, 5.0, // 25 + 5.0, 5.0, 5.0, 5.0, 5.0, // 30 + 5.0, 5.0, 5.0, 5.0, 5.0 // 35 +}; // in mm, one value for each plate, core and fin thickness +vector targetPS_plateLugThickness = { + 6.0, 6.0, 6.0, 6.0, 6.0, // 5 + 6.0, 6.0, 6.0, 6.0, 6.0, // 10 + 6.0, 6.0, 6.0, 6.0, 6.0, // 15 + 6.0, 6.0, 6.0, 6.0, 6.0, // 20 + 6.0, 6.0, 6.0, 6.0, 6.0, // 25 + 6.0, 6.0, 6.0, 6.0, 6.0, // 30 + 6.0, 6.0, 6.0, 6.0, 6.0 // 35 +}; // in mm, one value for each plate, note lugs are flush with the fins on the proton downstream side +bool targetPS_addFilletToPlateCore = true; // whether to add fillet to the intersection between plate core and fins +bool targetPS_addFilletToPlateLug = true; // whether to add fillet to the intersection between fins and lugs +// Fillets between fins and lugs are less important than the fillets between the core and fins, since the lugs are not in the direct path of the beam. +// Fillets due to lug thickness larger than the plate thickness is not added, due to the complexity of the geometry and the small effect on the physics. +double targetPS_plateFilletRadius = 1.0; // mm + +// Rods holding plates, which are attached to the end rings. +string targetPS_rodMaterial = "Inconel718"; +double targetPS_rodRadius = 1.5; // mm +// Rod length will be the sum of lug thicknesses and the spacer thicknesses, to be determined in the code. +// Actual rod is longer than this since it inserts into the end rings, but for the geometry reconstruction, that part will be treated as part of the end ring. + +// Spacer geometry, which is the same for all spacers. Added to both upstream and downstream end of the rods. +string targetPS_spacerMaterial = "Inconel718"; +double targetPS_spacerHalfLength = 1.5; // mm, full length 3.0 mm +double targetPS_spacerOuterRadius = 3.0; // mm +double targetPS_spacerInnerRadius = 1.55; // mm + +// End ring / support ring +string targetPS_supportRingMaterial = "Inconel718"; +double targetPS_supportRingLength = 8.1; // full length in mm +double targetPS_supportRingInnerRadius = 15.0; // mm +double targetPS_supportRingOuterRadius = 17.0; // mm +double targetPS_supportRingLugOuterRadius = 3.0; // mm +// Lug distance to center is the same as the plate lug distance to center, not defined as a separate parameter. +// Lug angles are constrained by the plate fin angles. +// End rings have chamfers at the sockets where the rods inserts (aligned with the lugs). The sockets, chamfers, and fillets at the bottom of the sockets are not modeled. +bool targetPS_addFilletToSupportRingLug = true; +double targetPS_supportRingLugFilletRadius = 1.0; +bool targetPS_addCutoutToSupportRing = true; // the six cylindrical cutouts at an angle +double targetPS_supportRingCutoutOffset = 4.78; // mm, distance from the center of the cutout to the side of the end ring closer to the target center. + // This value is a must-have. If no cutout is added, this determines where the support wire is attached from the inner edge of the end ring. +int targetPS_nSupportRingCutouts = 3; +vector targetPS_supportRingCutoutAngles = {15., 135., 255.}; // degrees + //{-45., 75., 195.}; these have spokes inside so not cutting out +double targetPS_supportRingCutoutInnerRadius = 2.06; // mm +double targetPS_supportRingCutoutTilt = 20.0; //degrees, angle of the cutout w.r.t the radial direction + +// Bicycle wheel parameters -- copied from the Hayman target +bool targetPS.supports.build = true; +double targetPS.supports.wheel.halfLength = 9.525; +double targetPS.supports.wheel.rOut = 196.85; +double targetPS.supports.wheel.rIn = 177.8; +string targetPS.supports.wheel.material = "G4_Al"; +int targetPS.supports.nSpokes = 3; +// Features on the wheel, changed according to docdb-30814 +vector targetPS.supports.features.angles = {-37., 83., 203.}; //degrees, center of the features +vector targetPS.supports.features.arcs = { 28., 28., 28.}; //degrees +vector targetPS.supports.features.rIns = {158.75, 158.75, 158.75}; //mm +//support rods in the wheel that the wires connect to +vector targetPS.supports.rods.angles = {-45., 75., 195.0}; //degrees +vector targetPS.supports.rods.halfLength = {70., 70., 70.}; //mm +vector targetPS.supports.rods.offset = {-35., -14.80, 38.62}; //mm, changed according to F10122537_A_DWG1 + //rod center to the channel centers +vector targetPS.supports.rods.pinOffset = {3.81, 3.81, 3.81}; //mm, according to F10121763_A_DWG1, pinholes used to fix + //the rods are off-center, need to add pinOffset/cos(rotY) to offset above +vector targetPS.supports.rods.radius = {9., 9., 9.}; //mm +vector targetPS.supports.rods.radialOffset = {177.8, 177.8, 177.8}; //mm +vector targetPS.supports.rods.wireOffset.downstream = {13., 13., 13.}; //mm, on the surface +vector targetPS.supports.rods.wireOffset.upstream = {13., 13., 13.}; //mm, calculated according to F10122537_A_DWG1 +//spokes connecting the target to the wheel rods +vector targetPS.supports.spokes.targetAngles.downstream = {-45., 75., 195.0}; //degrees +vector targetPS.supports.spokes.targetAngles.upstream = {-45., 75., 195.0}; //degrees +double targetPS.supports.spokes.diameter = 2.; +string targetPS.supports.spokes.material = "Inconel718"; + +// +// This tells emacs to view this file in c++ mode. +// Local Variables: +// mode:c++ +// End: diff --git a/Mu2eG4/geom/geom_common.txt b/Mu2eG4/geom/geom_common.txt index 648a3268fc..b9c624ecf9 100644 --- a/Mu2eG4/geom/geom_common.txt +++ b/Mu2eG4/geom/geom_common.txt @@ -3,7 +3,7 @@ // The file that is included will be time dependent. // -#include "Offline/Mu2eG4/geom/geom_run1_a.txt" +#include "Offline/Mu2eG4/geom/geom_run1_a_stickman.txt" // This tells emacs to view this file in c++ mode. // Local Variables: diff --git a/Mu2eG4/geom/geom_run1_a_stickman.txt b/Mu2eG4/geom/geom_run1_a_stickman.txt new file mode 100644 index 0000000000..47ed2d2bc8 --- /dev/null +++ b/Mu2eG4/geom/geom_run1_a_stickman.txt @@ -0,0 +1,118 @@ +// Based off Mu2eG4/geom/geom_run1_a.txt +// Updates made to include the new Stickman muon production target + +string detector.name = "g4geom_v00"; + +bool hasHall = true; +bool hasTarget = true; +bool hasProtonAbsorber = true; +bool hasTSdA = true; +bool hasExternalShielding = true; +bool hasDiskCalorimeter = true; +bool hasBeamline = true; +bool hasVirtualDetector = true; // some components, e.g. ProtonAbsorber assume vd presence now; +bool hasCosmicRayShield = true; +bool hasSTM = true; +bool hasMBS = true; // note the two subcomponents, see mbs section below; + // no MBS implies no downstream hole in Cosmic Ray Passive Shield + // Magnetic field may be affected as well +bool hasPTM = true; + + + + +#include "Offline/Mu2eG4/geom/g4_visOptions.txt" + +//------------------------------------------- +// Mu2e geometry includes +//------------------------------------------- + +// X-offset of the PS(+x) and DS(-x) from the Mu2e origin. +// The origin of the detector coordinate system is on the DS axis at the specified z. +double mu2e.solenoidOffset = 3904.; // mm +double mu2e.detectorSystemZ0 = 10171.; // mm G4BL: (17730-7292=9801 mm) + +#include "Offline/Mu2eG4/geom/mu2eWorld.txt" +// mu2eHall.txt should be used with protonBeamDump_v02.txt, below +//#include "Mu2eG4/geom/mu2eHall.txt" +// whereas mu2eHall_v*.txt should be used with protonBeamDump_v03.txt, below +#include "Offline/Mu2eG4/geom/mu2eHall_v04.txt" + +// Solenoids +#include "Offline/Mu2eG4/geom/DetectorSolenoid_v05.txt" +#include "Offline/Mu2eG4/geom/DSShielding_v03.txt" +#include "Offline/Mu2eG4/geom/ProductionSolenoid_v02.txt" +#include "Offline/Mu2eG4/geom/psEnclosure_v05.txt" +#include "Offline/Mu2eG4/geom/PSShield_v06.txt" +#include "Offline/Mu2eG4/geom/PSExternalShielding_v01.txt" +#include "Offline/Mu2eG4/geom/TransportSolenoid_v08.txt" + +// External Shielding +#include "Offline/Mu2eG4/geom/reduced_ExtShieldUpstream_v07.txt" +#include "Offline/Mu2eG4/geom/reduced_ExtShieldDownstream_v06.txt" +#include "Offline/Mu2eG4/geom/Saddle_v03.txt" +#include "Offline/Mu2eG4/geom/Pipe_v04.txt" +#include "Offline/Mu2eG4/geom/ElectronicRack_v02.txt" + +// #include "Mu2eG4/geom/stoppingTargetHoles_DOE_review_2017.txt" // 37 foil muon stopping target with holes +#include "Offline/Mu2eG4/geom/stoppingTargetHoles_v02.txt" // 37 foil muon stopping target with holes and thicker foils (0.1 mm -> 0.1056 mm) + +#include "Offline/Mu2eG4/geom/TSdA_v02.txt" +#include "Offline/Mu2eG4/geom/muonBeamStop_v08.txt" + +#include "Offline/Mu2eG4/geom/STM_v09.txt" // (muon) stopping target monitor + +// Proton Absorber +#include "Offline/Mu2eG4/geom/protonAbsorber_cylindrical_v04.txt" +#include "Offline/Mu2eG4/geom/degrader_v02.txt" // pion degrader. Off by default + +#include "Offline/Mu2eG4/geom/ProductionTarget_Stickman_v1_0.txt" +#include "Offline/Mu2eG4/geom/protonBeamDump_v03.txt" +#include "Offline/Mu2eG4/geom/extmon_fnal_v02.txt" + +// as-built panel numbering +#include "Offline/Mu2eG4/geom/tracker_v7.txt" + +// Crystal calorimeter +#include "Offline/Mu2eG4/geom/calorimeter_CsI_v2.txt" + +//CRV counters +#include "Offline/Mu2eG4/geom/crv_counters_run1a_v01.txt" + +// Production target beam-scanning detectors +#include "Offline/Mu2eG4/geom/PTM_v02.txt" + +//--------------------------------------- +// Virtual detectors +//--------------------------------------- +double vd.halfLength = 0.01; //mm +int vd.verbosityLevel = 0; +bool vd.visible = true; +bool vd.solid = false; + +// // VD right in front of a hall wall +// double vd.ExtMonCommonPlane.z = -11999.99; + + +//--------------------------------------- +// Region visualization +//--------------------------------------- +#include "Offline/Mu2eG4/geom/visualization_regions.txt" + + +// patch to change East DS hatch block height +double building.dsArea.hatchblock.E.offsetFromFloorSurface.y = 6400.80; +double building.dsArea.hatchblock.E.yHalfThickness = 457.20; + +// +// +// End notes: +// +// 1) Sources of information: +// +// +// +// This tells emacs to view this file in c++ mode. +// Local Variables: +// mode:c++ +// End: diff --git a/Mu2eG4/inc/constructTargetPS.hh b/Mu2eG4/inc/constructTargetPS.hh index 3e9c60ca78..cbc7ef159c 100644 --- a/Mu2eG4/inc/constructTargetPS.hh +++ b/Mu2eG4/inc/constructTargetPS.hh @@ -18,6 +18,11 @@ namespace mu2e { void constructTargetPS(VolumeInfo const & parent, SimpleConfig const & _config); + // Target-type-specific construction functions + void constructTier1Target(VolumeInfo const & parent, SimpleConfig const & _config); + void constructHaymanTarget(VolumeInfo const & parent, SimpleConfig const & _config); + void constructStickmanTarget(VolumeInfo const & parent, SimpleConfig const & _config); + } #endif /* Mu2eG4_constructTargetPS_hh */ diff --git a/Mu2eG4/src/constructPS.cc b/Mu2eG4/src/constructPS.cc index 22d0bdf625..c5955ecd11 100644 --- a/Mu2eG4/src/constructPS.cc +++ b/Mu2eG4/src/constructPS.cc @@ -9,6 +9,7 @@ // C++ includes #include +#include // Mu2e includes. #include "Offline/BeamlineGeom/inc/Beamline.hh" @@ -356,22 +357,25 @@ namespace mu2e { } // std::cout << "inside " << __func__ << psVacuumParams.originInMu2e() << " " << _hallOriginInMu2e << std::endl; - // Build the production target. - - if ( targetPS_model == "MDC2018" ){ - verbosityLevel> 0 && std::cout << __func__ << "MDC 2018 target" << std::endl; - constructTargetPS(psVacuumInfo, _config ); - } else if ( targetPS_model == "HaymanLowerDensity" ){ - verbosityLevel> 0 && std::cout << __func__ << "HaymanLowerDensity target" << std::endl; - constructHaymanRings(psVacuumInfo, _config); - } else if (targetPS_model == "Hayman_v_2_0"){ - verbosityLevel> 0 && std::cout << __func__ << "Hayman 2.0 target" << std::endl; - constructTargetPS(psVacuumInfo, _config ); - } else{ + // Build the production target using a dispatch registry to avoid growing if-else chains + using ConstructorFunction = void (*)(const VolumeInfo&, const SimpleConfig&); + static const std::map targetConstructors = { + {"MDC2018", &constructTargetPS}, + {"HaymanLowerDensity", &constructHaymanRings}, + {"Hayman_v_2_0", &constructTargetPS}, + {"Stickman_v_1_0", &constructTargetPS} + }; + + auto it = targetConstructors.find(targetPS_model); + if (it != targetConstructors.end()) { + if (verbosityLevel > 0) { + std::cout << __func__ << " constructing target model: " << targetPS_model << std::endl; + } + it->second(psVacuumInfo, _config); + } else { throw cet::exception("CONFIG") << "In constructPS.cc unrecognized production target model name: " - << targetPS_model - << "\n"; + << targetPS_model << "\n"; } // FIXME: make unconditional diff --git a/Mu2eG4/src/constructTargetPS.cc b/Mu2eG4/src/constructTargetPS.cc index a101584893..42e5625a3c 100644 --- a/Mu2eG4/src/constructTargetPS.cc +++ b/Mu2eG4/src/constructTargetPS.cc @@ -6,6 +6,7 @@ #include #include #include +#include // art includes #include "art/Framework/Services/Registry/ServiceDefinitionMacros.h" @@ -48,6 +49,7 @@ #include "Geant4/G4Tubs.hh" #include "Geant4/G4VPhysicalVolume.hh" #include "Geant4/G4PVPlacement.hh" +#include "Geant4/G4ExtrudedSolid.hh" #include "Geant4/G4UnionSolid.hh" #include "Geant4/G4IntersectionSolid.hh" @@ -77,525 +79,794 @@ namespace mu2e { finishNesting(boxWithTubsCutoutInfo,material,rotAngle,translation,parent.logical,copyNo,color,lookupToken,verbose); } + // Registry-based dispatcher for target model construction void constructTargetPS(VolumeInfo const & parent, SimpleConfig const & _config) { + // Registry mapping target versions to construction functions + using TargetConstructor = void (*)(VolumeInfo const &, SimpleConfig const &); + static const std::map targetRegistry = { + {ProductionTargetMaker::tier1, &constructTier1Target}, + {ProductionTargetMaker::hayman_v_2_0, &constructHaymanTarget}, + {ProductionTargetMaker::stickman_v_1_0, &constructStickmanTarget} + }; + + int targetVersion = _config.getInt("targetPS_version"); + auto it = targetRegistry.find(targetVersion); + + if (it != targetRegistry.end()) { + it->second(parent, _config); + } else { + throw cet::exception("CONFIG") + << "In constructTargetPS: unrecognized production target version = " + << targetVersion << "\n"; + } + } //end constructTargetPS + + // ========== Target-specific construction functions ========== + + void constructTier1Target(VolumeInfo const & parent, SimpleConfig const & _config) { Mu2eG4Helper & _helper = *(art::ServiceHandle()); AntiLeakRegistry & reg = _helper.antiLeakRegistry(); - // - // which target am I building? - if (_config.getInt("targetPS_version") == ProductionTargetMaker::tier1){ - int verbosityLevel = _config.getInt("PS.verbosityLevel"); + int verbosityLevel = _config.getInt("PS.verbosityLevel"); + + //G4Material* psVacVesselMaterial = findMaterialOrThrow(psVacVesselInnerParams.materialName()); + + verbosityLevel >0 && + cout << __func__ << " verbosityLevel : " << verbosityLevel << endl; + + //bool psVacuumSensitive = _config.getBool("PS.Vacuum.Sensitive", false); + + G4ThreeVector _hallOriginInMu2e = parent.centerInMu2e(); + + + // Build the production target. + GeomHandle tgt; + + double const _clamp_supWheel_rOut = _config.getDouble("clamp_supWheel_rOut"); + double const _clamp_supWheel_halfLength = _config.getDouble("clamp_supWheel_halfLength"); + double const _supWheel_trgtPS_halfLength = _config.getDouble("supWheel_trgtPS_halfLength"); + + double envHalfLength = 2.0*_clamp_supWheel_halfLength+_supWheel_trgtPS_halfLength; + if (tgt->envHalfLength()>envHalfLength) { envHalfLength = tgt->envHalfLength(); } + + + TubsParams prodTargetMotherParams( 0., _clamp_supWheel_rOut, envHalfLength); + VolumeInfo prodTargetMotherInfo = nestTubs( "ProductionTargetMother", + prodTargetMotherParams, + parent.logical->GetMaterial(), + 0, + tgt->position() - parent.centerInMu2e(), + parent, + 0, + G4Colour::Blue(), + "PS" + ); + if(verbosityLevel > 0) { + G4cout << "inside clause 1 of " << __func__ << G4endl; + G4cout << "target position and hall origin = " << tgt->position() << "\n" << _hallOriginInMu2e << " " << + parent.centerInMu2e() << G4endl; + G4cout << "supWheel and envHalfLength = " << _clamp_supWheel_rOut << "\n" << + envHalfLength << G4endl; + } + + double const _supWheel_trgtPS_rIn = _config.getDouble("supWheel_trgtPS_rIn"); + double const _supWheel_trgtPS_rOut = _config.getDouble("supWheel_trgtPS_rOut"); + G4Material* _suppWheelMaterial = findMaterialOrThrow(_config.getString("supWheel_trgtPS_material")); + + G4GeometryOptions* geomOptions = art::ServiceHandle()->geomOptions(); + geomOptions->loadEntry( _config, "ProductionTargetSupportWheel", "supWheel_trgtPS" ); + + G4ThreeVector _loclCenter(0.0,0.0,0.0); + + TubsParams suppWheelParams( _supWheel_trgtPS_rIn, _supWheel_trgtPS_rOut, _supWheel_trgtPS_halfLength); + VolumeInfo suppWheelInfo = nestTubs( "ProductionTargetSupportWheel", + suppWheelParams, + _suppWheelMaterial, + 0, + _loclCenter, + prodTargetMotherInfo, + 0, + G4Colour::Gray() + ); - //G4Material* psVacVesselMaterial = findMaterialOrThrow(psVacVesselInnerParams.materialName()); + double const _clamp_supWheel_rIn = _config.getDouble("clamp_supWheel_rIn"); + G4Material* _clampSpWheelMaterial = findMaterialOrThrow(_config.getString("clamp_supWheel_material")); + geomOptions->loadEntry( _config, "ClampSupportWheel_R", "clamp_supWheel" ); - verbosityLevel >0 && - cout << __func__ << " verbosityLevel : " << verbosityLevel << endl; + G4ThreeVector _clampPosR(0.0,0.0,_supWheel_trgtPS_halfLength+_clamp_supWheel_halfLength); - //bool psVacuumSensitive = _config.getBool("PS.Vacuum.Sensitive", false); + TubsParams clampSpWheelParams( _clamp_supWheel_rIn, _clamp_supWheel_rOut, _clamp_supWheel_halfLength); + VolumeInfo clampSpWheelInfoR = nestTubs( "ClampSupportWheel_R", + clampSpWheelParams, + _clampSpWheelMaterial, + 0, + _clampPosR, + prodTargetMotherInfo, + 0, + G4Colour::Gray() + ); - G4ThreeVector _hallOriginInMu2e = parent.centerInMu2e(); + // G4ThreeVector _clampPosL(0.0,0.0,-(_supWheel_trgtPS_halfLength+_clamp_supWheel_halfLength)); + // + // geomOptions->loadEntry( _config, "ClampSupportWheel_L", "clamp_supWheel" ); + // VolumeInfo clampSpWheelInfoL = nestTubs( "ClampSupportWheel_L", + // clampSpWheelParams, + // _clampSpWheelMaterial, + // 0, + // _clampPosL, + // prodTargetMotherInfo, + // 0, + // G4Colour::Gray() + // ); + + TubsParams prodTargetParams( 0., tgt->rOut(), tgt->halfLength()); + + G4Material* prodTargetMaterial = findMaterialOrThrow(_config.getString("targetPS_materialName")); + + geomOptions->loadEntry( _config, "ProductionTarget", "targetPS" ); + + VolumeInfo prodTargetInfo = nestTubs( "ProductionTarget", + prodTargetParams, + prodTargetMaterial, + &tgt->productionTargetRotation(), + _loclCenter, + prodTargetMotherInfo, + 0, + G4Colour::Magenta() + ); + if(verbosityLevel > 0) + G4cout << "_loclCenter for Tier 1 target = " << _loclCenter << G4endl; + // Now add fins for version 2 + if ( tgt->version() > 1 ) { + // Length of fins must be calculated: + + double finHalfLength = (tgt->halfLength() * 2.0 - tgt->hubDistUS() + - tgt->hubDistDS() - tgt->hubLenUS() - tgt->hubLenDS())/2.0; // on the inner side, + // adjacent to the main target body. + + // Use the steeper of the two hub angles to angle the ends of the fin + double theAngle = tgt->hubAngleUS(); + if ( tgt->hubAngleDS() > theAngle ) theAngle = tgt->hubAngleDS(); + + double finHalfLengthOut = finHalfLength + tgt->finHeight()*std::cos(theAngle*CLHEP::degree); + + G4Trd * myTrd = new G4Trd("FinTrapezoid", + finHalfLength, finHalfLengthOut, + tgt->finThickness()/2.0, tgt->finThickness()/2.0, + tgt->finHeight()/2.0); + + // std::vector finDims = {tgt->finThickness()/2.0,tgt->finHeight()/2.0,finHalfLength}; + double finZoff = (tgt->hubDistUS() - tgt->hubDistDS() + tgt->hubLenUS() - tgt->hubLenDS())/2.0; // z-offset for fin + double rToFin = tgt->rOut()+tgt->finHeight()/2.0+0.1; + + double xMove = finZoff * std::sin(tgt->productionTargetRotation().theta()); + CLHEP::Hep3Vector finOffset1(xMove,-rToFin,finZoff); + CLHEP::Hep3Vector finOffset2(rToFin*std::cos(M_PI/6.0)+xMove,rToFin*std::sin(M_PI/6.0),finZoff-1.0); + CLHEP::Hep3Vector finOffset3(-rToFin*std::cos(M_PI/6.0)+xMove-0.1,rToFin*std::sin(M_PI/6.0)+0.02,finZoff+0.2); + CLHEP::HepRotation* rotFinBase = reg.add(CLHEP::HepRotation(CLHEP::HepRotation::IDENTITY)); + rotFinBase->rotateX(90.0*CLHEP::degree); + rotFinBase->rotateZ(90.0*CLHEP::degree); + + CLHEP::HepRotation* rotFin1 = reg.add(CLHEP::HepRotation((*rotFinBase)*tgt->productionTargetRotation())); + CLHEP::HepRotation* rotFin2 = reg.add(CLHEP::HepRotation((*rotFinBase)*tgt->productionTargetRotation())); + CLHEP::HepRotation* rotFin3 = reg.add(CLHEP::HepRotation((*rotFinBase)*tgt->productionTargetRotation())); + rotFin1->rotateX(M_PI); + rotFin2->rotateX(M_PI/3.0); + rotFin3->rotateX(5.0*M_PI/3.0); + VolumeInfo fin1Vol("ProductionTargetFin1", + _loclCenter, + prodTargetMotherInfo.centerInWorld); + + VolumeInfo fin2Vol("ProductionTargetFin2", + _loclCenter, + prodTargetMotherInfo.centerInWorld); + + VolumeInfo fin3Vol("ProductionTargetFin3", + _loclCenter, + prodTargetMotherInfo.centerInWorld); + + fin1Vol.solid = myTrd; + fin2Vol.solid = myTrd; + fin3Vol.solid = myTrd; + + finishNesting( fin1Vol, + prodTargetMaterial, + rotFin1, + finOffset1, + prodTargetMotherInfo.logical, + 0, + G4Colour::Magenta(), + "PS" + ); + finishNesting( fin2Vol, + prodTargetMaterial, + rotFin2, + finOffset2, + prodTargetMotherInfo.logical, + 0, + G4Colour::Magenta(), + "PS" + ); + + finishNesting( fin3Vol, + prodTargetMaterial, + rotFin3, + finOffset3, + prodTargetMotherInfo.logical, + 0, + G4Colour::Magenta(), + "PS" + ); - // Build the production target. - GeomHandle tgt; + } - double const _clamp_supWheel_rOut = _config.getDouble("clamp_supWheel_rOut"); - double const _clamp_supWheel_halfLength = _config.getDouble("clamp_supWheel_halfLength"); - double const _supWheel_trgtPS_halfLength = _config.getDouble("supWheel_trgtPS_halfLength"); - double envHalfLength = 2.0*_clamp_supWheel_halfLength+_supWheel_trgtPS_halfLength; - if (tgt->envHalfLength()>envHalfLength) { envHalfLength = tgt->envHalfLength(); } + // Using the old terms "right" and "left" to mean "downstream" (DS) + // and "upstream" (US), respectively. - TubsParams prodTargetMotherParams( 0., _clamp_supWheel_rOut, envHalfLength); - VolumeInfo prodTargetMotherInfo = nestTubs( "ProductionTargetMother", - prodTargetMotherParams, - parent.logical->GetMaterial(), + Polycone const & pHubRgtParams = *tgt->getHubsRgtPtr(); + VolumeInfo prodTargetHubRgtInfo = nestPolycone("ProductionTargetHubRgt", + pHubRgtParams.getPolyconsParams(), + findMaterialOrThrow(pHubRgtParams.materialName()), + &tgt->productionTargetRotation(), + pHubRgtParams.originInMu2e(), + prodTargetMotherInfo, 0, - tgt->position() - parent.centerInMu2e(), - parent, + G4Colour::Magenta(), + "ProductionTarget" + ); + + Polycone const & pHubLftParams = *tgt->getHubsLftPtr(); + VolumeInfo prodTargetHubLftInfo = nestPolycone("ProductionTargetHubLft", + pHubLftParams.getPolyconsParams(), + findMaterialOrThrow(pHubLftParams.materialName()), + &tgt->productionTargetRotation(), + pHubLftParams.originInMu2e(), + prodTargetMotherInfo, 0, - G4Colour::Blue(), - "PS" + G4Colour::Magenta(), + "ProductionTarget" ); - if(verbosityLevel > 0) { - G4cout << "inside clause 1 of " << __func__ << G4endl; - G4cout << "target position and hall origin = " << tgt->position() << "\n" << _hallOriginInMu2e << " " << - parent.centerInMu2e() << G4endl; - G4cout << "supWheel and envHalfLength = " << _clamp_supWheel_rOut << "\n" << - envHalfLength << G4endl; - } - double const _supWheel_trgtPS_rIn = _config.getDouble("supWheel_trgtPS_rIn"); - double const _supWheel_trgtPS_rOut = _config.getDouble("supWheel_trgtPS_rOut"); - G4Material* _suppWheelMaterial = findMaterialOrThrow(_config.getString("supWheel_trgtPS_material")); - - G4GeometryOptions* geomOptions = art::ServiceHandle()->geomOptions(); - geomOptions->loadEntry( _config, "ProductionTargetSupportWheel", "supWheel_trgtPS" ); - - G4ThreeVector _loclCenter(0.0,0.0,0.0); - - TubsParams suppWheelParams( _supWheel_trgtPS_rIn, _supWheel_trgtPS_rOut, _supWheel_trgtPS_halfLength); - VolumeInfo suppWheelInfo = nestTubs( "ProductionTargetSupportWheel", - suppWheelParams, - _suppWheelMaterial, - 0, - _loclCenter, - prodTargetMotherInfo, - 0, - G4Colour::Gray() - ); - - double const _clamp_supWheel_rIn = _config.getDouble("clamp_supWheel_rIn"); - G4Material* _clampSpWheelMaterial = findMaterialOrThrow(_config.getString("clamp_supWheel_material")); - geomOptions->loadEntry( _config, "ClampSupportWheel_R", "clamp_supWheel" ); - - G4ThreeVector _clampPosR(0.0,0.0,_supWheel_trgtPS_halfLength+_clamp_supWheel_halfLength); - - TubsParams clampSpWheelParams( _clamp_supWheel_rIn, _clamp_supWheel_rOut, _clamp_supWheel_halfLength); - VolumeInfo clampSpWheelInfoR = nestTubs( "ClampSupportWheel_R", - clampSpWheelParams, - _clampSpWheelMaterial, - 0, - _clampPosR, - prodTargetMotherInfo, - 0, - G4Colour::Gray() - ); - - // G4ThreeVector _clampPosL(0.0,0.0,-(_supWheel_trgtPS_halfLength+_clamp_supWheel_halfLength)); + CLHEP::Hep3Vector zax(0,0,1); + double spokeRad = 0.5*_config.getDouble("targetPS_Spoke_diameter"); + G4Material* spokeMaterial = findMaterialOrThrow(_config.getString("targetPS_Spoke_materialName")); + double hub_overhang_angle = _config.getDouble("targetPS_Hub_overhang_angle")*CLHEP::deg; + double normalOverhangRgt = CLHEP::halfpi+hub_overhang_angle; + double normalOverhangLft = CLHEP::halfpi-hub_overhang_angle; + CLHEP::HepRotation invTrgtRot = tgt->productionTargetRotation().inverse(); + + std::map const & anchoringPntsRgt = tgt->anchoringPntsRgt(); + + int iSpk(0); + + // Create variable to avoid multiple look-up + const bool prodTargetVisible = geomOptions->isVisible( "ProductionTarget" ); + const bool prodTargetSolid = geomOptions->isSolid ( "ProductionTarget" ); + const bool forceAuxEdgeVisible = geomOptions->forceAuxEdgeVisible( "ProductionTarget" ); + const bool placePV = geomOptions->placePV( "ProductionTarget" ); + const bool doSurfaceCheck = geomOptions->doSurfaceCheck( "ProductionTarget" ); + + for (std::map::const_iterator iPnt=anchoringPntsRgt.begin(); iPnt!=anchoringPntsRgt.end(); ++iPnt) { + double tmpAngle = iPnt->first * CLHEP::deg; + CLHEP::Hep3Vector normalax(std::cos(tmpAngle),std::sin(tmpAngle),0.0); + CLHEP::Hep3Vector tmpAnchPnt = _supWheel_trgtPS_rIn*normalax; + CLHEP::Hep3Vector tmpDirVec = tmpAnchPnt-iPnt->second; + CLHEP::Hep3Vector tmpMidPnt(0.5*(tmpAnchPnt.x()+iPnt->second.x()),0.5*(tmpAnchPnt.y()+iPnt->second.y()),0.5*(tmpAnchPnt.z()+iPnt->second.z())); + double tmpSpokeLength = tmpDirVec.mag(); + tmpDirVec = tmpDirVec.unit(); + CLHEP::Hep3Vector tmpRotAxis = tmpDirVec.cross(zax); + double rotAngle = tmpDirVec.angle(zax); + double deepAngleWheel = tmpDirVec.angle(normalax); + + CLHEP::Hep3Vector normalaxHub(std::sin(normalOverhangRgt),0,std::cos(normalOverhangRgt)); + normalaxHub.rotateZ(tmpAngle); + normalaxHub.transform(invTrgtRot); + double deepAngleHub = tmpDirVec.angle(normalaxHub); + + double fixOverlapWheel = spokeRad*std::tan(deepAngleWheel); + double fixOverlapHub = spokeRad*std::tan(deepAngleHub); + tmpSpokeLength-=(fixOverlapWheel+fixOverlapHub); + tmpMidPnt += 0.5*(fixOverlapHub-fixOverlapWheel)*tmpDirVec; + + TubsParams iSpokeParams( 0.0, spokeRad, 0.5*tmpSpokeLength); + std::stringstream iSpokeName; + iSpokeName<<"ProductionTargetSpokeRgt_"< const & anchoringPntsLft = tgt->anchoringPntsLft(); + + iSpk=0; + for (std::map::const_iterator iPnt=anchoringPntsLft.begin(); iPnt!=anchoringPntsLft.end(); ++iPnt) { + double tmpAngle = iPnt->first * CLHEP::deg; + CLHEP::Hep3Vector normalax(std::cos(tmpAngle),std::sin(tmpAngle),0.0); + CLHEP::Hep3Vector tmpAnchPnt = _supWheel_trgtPS_rIn*normalax; + CLHEP::Hep3Vector tmpDirVec = tmpAnchPnt-iPnt->second; + CLHEP::Hep3Vector tmpMidPnt(0.5*(tmpAnchPnt.x()+iPnt->second.x()),0.5*(tmpAnchPnt.y()+iPnt->second.y()),0.5*(tmpAnchPnt.z()+iPnt->second.z())); + double tmpSpokeLength = tmpDirVec.mag(); + tmpDirVec = tmpDirVec.unit(); + CLHEP::Hep3Vector tmpRotAxis = tmpDirVec.cross(zax); + double rotAngle = tmpDirVec.angle(zax); + double deepAngleWheel = tmpDirVec.angle(normalax); + + CLHEP::Hep3Vector normalaxHub(std::sin(normalOverhangLft),0,std::cos(normalOverhangLft)); + normalaxHub.rotateZ(tmpAngle); + normalaxHub.transform(invTrgtRot); + double deepAngleHub = tmpDirVec.angle(normalaxHub); + //CLHEP::HepRotation *tmpSpokeRot = reg.add(CLHEP::HepRotation(tmpRotAxis,rotAngle)); + + double fixOverlapWheel = spokeRad*std::tan(deepAngleWheel); + double fixOverlapHub = spokeRad*std::tan(deepAngleHub); + tmpSpokeLength-=(fixOverlapWheel+fixOverlapHub); + tmpMidPnt += 0.5*(fixOverlapHub-fixOverlapWheel)*tmpDirVec; + + TubsParams iSpokeParams( 0.0, spokeRad, 0.5*tmpSpokeLength); + std::stringstream iSpokeName; + iSpokeName<<"ProductionTargetSpokeLft_"<()); + AntiLeakRegistry & reg = _helper.antiLeakRegistry(); + + int verbosityLevel = _config.getInt("PSHayman.verbosityLevel"); + verbosityLevel >0 && + cout << __func__ << " verbosityLevel on Hayman2.0 : " << verbosityLevel << endl; + G4ThreeVector _hallOriginInMu2e = parent.centerInMu2e(); + // Create variable to avoid multiple look-up + + G4GeometryOptions* geomOptions = art::ServiceHandle()->geomOptions(); + + bool prodTargetVisible = geomOptions->isVisible( "ProductionTarget" ); + bool prodTargetSolid = geomOptions->isSolid ( "ProductionTarget" ); + bool forceAuxEdgeVisible = geomOptions->forceAuxEdgeVisible( "ProductionTarget" ); + bool placePV = geomOptions->placePV( "ProductionTarget" ); + bool doSurfaceCheck = geomOptions->doSurfaceCheck( "ProductionTarget" ); + // + + // begin all names with ProductionTarget so when we build sensitive detectors we can make them all + // sensitive at once with LVname.find("ProductionTarget") !=std::string::npos + // + // Build the production target. + GeomHandle tgt; + G4Material* prodTargetCoreMaterial = findMaterialOrThrow(tgt->targetCoreMaterial()); + G4Material* prodTargetFinMaterial = findMaterialOrThrow(tgt->targetFinMaterial()); + G4Material* prodTargetSupportRingMaterial = findMaterialOrThrow(tgt->supportRingMaterial()); + + TubsParams prodTargetMotherParams( 0. + ,tgt->productionTargetMotherOuterRadius() + ,tgt->productionTargetMotherHalfLength()); + + G4ThreeVector _loclCenter(0.0,0.0,0.0); + G4RotationMatrix* targetRotation = reg.add(G4RotationMatrix(tgt->productionTargetRotation().inverse())); + if (verbosityLevel > 2){G4cout << __PRETTY_FUNCTION__ << "target rotation = " << *targetRotation << G4endl;} + VolumeInfo prodTargetMotherInfo = nestTubs( "ProductionTargetMother", + prodTargetMotherParams, + parent.logical->GetMaterial(), + 0, + tgt->haymanProdTargetPosition() - parent.centerInMu2e(), + parent, + 0, + G4Colour::Blue(), + "PS" + ); + + if (verbosityLevel > 0){ + G4cout << "target position and hall origin = " + << tgt->haymanProdTargetPosition() << "\n" << + _hallOriginInMu2e << " " << parent.centerInMu2e() << G4endl; + } + if (verbosityLevel > 2){ + G4cout << __PRETTY_FUNCTION__ << "created prodTargetMotherInfo " + << tgt->productionTargetMotherOuterRadius() + << " " <productionTargetMotherHalfLength() << G4endl; + } + // + // construct each section + + int numberOfSections = tgt->numberOfTargetSections(); + + + + // + // first build the core + + // + // start at most negative end (technically could start at center but this is simpler since target is not symmetric about its center + // this variable gets incremented every time we add anything + double _currentZ = _loclCenter.z() - tgt->halfHaymanLength(); + if (verbosityLevel > 2){G4cout << __PRETTY_FUNCTION__ << " current Z starts at " << _currentZ << G4endl;} + // + // running z locations for placing elements + CLHEP::Hep3Vector _segmentCenter(0.,0.,0.); + CLHEP::Hep3Vector _startingSegmentCenter(0.,0.,0.); + double targetRadius = tgt->rOut(); + // + // keeping track of these for sensitive volume creation + int coreCopyNumber = 0; + int finCopyNumber = 0; + int finTopCopyNumber = 0; + + double angularSize = std::asin((tgt->haymanFinThickness()/2.)/targetRadius); + double dSphi = M_PI/2. - angularSize; + double dPphi = 2.*angularSize; + + // + // need this for a) drawing subtraction volumes and b) rectangular cutouts out of circular support rings, etc. otherwise visualizer gets + // confused and how can G4 know which of two volumes some plane, for example, belongs to? + constexpr double magicOffset = 0.0001; + + // + //other constants + double cutoutHeightAboveTarget = tgt->supportRingInnerRadius() + tgt->supportRingCutoutThickness()/2. - magicOffset; + + CLHEP::HepRotation* rotRing = reg.add(CLHEP::HepRotation(tgt->productionTargetRotation())); + + // + // we're going to use these over and over; save time, memory, make cleaner and put in registry. save typing with currentFinAngles + // fins are constructed at 90^o relative to boxes and it's just cleaner to do it this way + std::vector rotFinsG4; + std::vector rotBoxesG4; + std::vector currentFinAngles; + for (int ithFin = 0; ithFin < tgt->nHaymanFins(); ++ithFin){ + rotFinsG4.emplace_back(reg.add(CLHEP::HepRotation(tgt->productionTargetRotation()))); + rotBoxesG4.emplace_back(reg.add(CLHEP::HepRotation::IDENTITY)); // - // geomOptions->loadEntry( _config, "ClampSupportWheel_L", "clamp_supWheel" ); - // VolumeInfo clampSpWheelInfoL = nestTubs( "ClampSupportWheel_L", - // clampSpWheelParams, - // _clampSpWheelMaterial, - // 0, - // _clampPosL, - // prodTargetMotherInfo, - // 0, - // G4Colour::Gray() - // ); - - TubsParams prodTargetParams( 0., tgt->rOut(), tgt->halfLength()); - - G4Material* prodTargetMaterial = findMaterialOrThrow(_config.getString("targetPS_materialName")); - - geomOptions->loadEntry( _config, "ProductionTarget", "targetPS" ); - - VolumeInfo prodTargetInfo = nestTubs( "ProductionTarget", - prodTargetParams, - prodTargetMaterial, - &tgt->productionTargetRotation(), - _loclCenter, - prodTargetMotherInfo, - 0, - G4Colour::Magenta() - ); - if(verbosityLevel > 0) - G4cout << "_loclCenter for Tier 1 target = " << _loclCenter << G4endl; - // Now add fins for version 2 - if ( tgt->version() > 1 ) { - // Length of fins must be calculated: - - double finHalfLength = (tgt->halfLength() * 2.0 - tgt->hubDistUS() - - tgt->hubDistDS() - tgt->hubLenUS() - tgt->hubLenDS())/2.0; // on the inner side, - // adjacent to the main target body. - - // Use the steeper of the two hub angles to angle the ends of the fin - double theAngle = tgt->hubAngleUS(); - if ( tgt->hubAngleDS() > theAngle ) theAngle = tgt->hubAngleDS(); - - double finHalfLengthOut = finHalfLength + tgt->finHeight()*std::cos(theAngle*CLHEP::degree); - - G4Trd * myTrd = new G4Trd("FinTrapezoid", - finHalfLength, finHalfLengthOut, - tgt->finThickness()/2.0, tgt->finThickness()/2.0, - tgt->finHeight()/2.0); - - // std::vector finDims = {tgt->finThickness()/2.0,tgt->finHeight()/2.0,finHalfLength}; - double finZoff = (tgt->hubDistUS() - tgt->hubDistDS() + tgt->hubLenUS() - tgt->hubLenDS())/2.0; // z-offset for fin - double rToFin = tgt->rOut()+tgt->finHeight()/2.0+0.1; - - double xMove = finZoff * sin(tgt->productionTargetRotation().theta()); - CLHEP::Hep3Vector finOffset1(xMove,-rToFin,finZoff); - CLHEP::Hep3Vector finOffset2(rToFin*cos(M_PI/6.0)+xMove,rToFin*sin(M_PI/6.0),finZoff-1.0); - CLHEP::Hep3Vector finOffset3(-rToFin*cos(M_PI/6.0)+xMove-0.1,rToFin*sin(M_PI/6.0)+0.02,finZoff+0.2); - CLHEP::HepRotation* rotFinBase = reg.add(CLHEP::HepRotation(CLHEP::HepRotation::IDENTITY)); - rotFinBase->rotateX(90.0*CLHEP::degree); - rotFinBase->rotateZ(90.0*CLHEP::degree); - - CLHEP::HepRotation* rotFin1 = reg.add(CLHEP::HepRotation((*rotFinBase)*tgt->productionTargetRotation())); - CLHEP::HepRotation* rotFin2 = reg.add(CLHEP::HepRotation((*rotFinBase)*tgt->productionTargetRotation())); - CLHEP::HepRotation* rotFin3 = reg.add(CLHEP::HepRotation((*rotFinBase)*tgt->productionTargetRotation())); - rotFin1->rotateX(M_PI); - rotFin2->rotateX(M_PI/3.0); - rotFin3->rotateX(5.0*M_PI/3.0); - VolumeInfo fin1Vol("ProductionTargetFin1", - _loclCenter, - prodTargetMotherInfo.centerInWorld); - - VolumeInfo fin2Vol("ProductionTargetFin2", - _loclCenter, - prodTargetMotherInfo.centerInWorld); - - VolumeInfo fin3Vol("ProductionTargetFin3", - _loclCenter, - prodTargetMotherInfo.centerInWorld); - - fin1Vol.solid = myTrd; - fin2Vol.solid = myTrd; - fin3Vol.solid = myTrd; - - finishNesting( fin1Vol, - prodTargetMaterial, - rotFin1, - finOffset1, - prodTargetMotherInfo.logical, - 0, - G4Colour::Magenta(), - "PS" - ); - - finishNesting( fin2Vol, - prodTargetMaterial, - rotFin2, - finOffset2, - prodTargetMotherInfo.logical, - 0, - G4Colour::Magenta(), - "PS" - ); - - finishNesting( fin3Vol, - prodTargetMaterial, - rotFin3, - finOffset3, - prodTargetMotherInfo.logical, - 0, - G4Colour::Magenta(), - "PS" - ); + // here I have a little confusing fix. The fin is built along the y-axis. But the rotation is given wrt the x axis. Hence I need to + // subtract off that 90^o . Then the - sign is CLHEP vs G4 + currentFinAngles.emplace_back(tgt->finAngles().at(ithFin)); + rotFinsG4.at(ithFin)->rotateZ(-(-M_PI/2. + currentFinAngles.at(ithFin))); + rotBoxesG4.at(ithFin)->rotateZ(-(-M_PI/2. + currentFinAngles.at(ithFin))); + if (verbosityLevel > 1){G4cout << "rotFin " << ithFin << " " << *rotFinsG4.at(ithFin) << G4endl;} + } - } + G4Box* cutoutBox = reg.add(G4Box("cutoutBox",tgt->supportRingCutoutThickness()/2.+ magicOffset, + tgt->haymanFinThickness() + magicOffset,tgt->supportRingCutoutLength()/2.+ magicOffset)); + // get real + for (int ithSection = 0; ithSection < numberOfSections; ++ithSection){ + int numberOfSegments = tgt->numberOfSegmentsPerSection().at(ithSection); + std::string startName = "ProductionTargetStartingCoreSection_" + std::to_string(ithSection+1); + // + // first place starting segment - // Using the old terms "right" and "left" to mean "downstream" (DS) - // and "upstream" (US), respectively. - - Polycone const & pHubRgtParams = *tgt->getHubsRgtPtr(); - VolumeInfo prodTargetHubRgtInfo = nestPolycone("ProductionTargetHubRgt", - pHubRgtParams.getPolyconsParams(), - findMaterialOrThrow(pHubRgtParams.materialName()), - &tgt->productionTargetRotation(), - pHubRgtParams.originInMu2e(), - prodTargetMotherInfo, - 0, - G4Colour::Magenta(), - "ProductionTarget" - ); - - Polycone const & pHubLftParams = *tgt->getHubsLftPtr(); - VolumeInfo prodTargetHubLftInfo = nestPolycone("ProductionTargetHubLft", - pHubLftParams.getPolyconsParams(), - findMaterialOrThrow(pHubLftParams.materialName()), - &tgt->productionTargetRotation(), - pHubLftParams.originInMu2e(), - prodTargetMotherInfo, - 0, - G4Colour::Magenta(), - "ProductionTarget" - ); - - CLHEP::Hep3Vector zax(0,0,1); - double spokeRad = 0.5*_config.getDouble("targetPS_Spoke_diameter"); - G4Material* spokeMaterial = findMaterialOrThrow(_config.getString("targetPS_Spoke_materialName")); - double hub_overhang_angle = _config.getDouble("targetPS_Hub_overhang_angle")*CLHEP::deg; - double normalOverhangRgt = CLHEP::halfpi+hub_overhang_angle; - double normalOverhangLft = CLHEP::halfpi-hub_overhang_angle; - CLHEP::HepRotation invTrgtRot = tgt->productionTargetRotation().inverse(); - - std::map const & anchoringPntsRgt = tgt->anchoringPntsRgt(); - - int iSpk(0); - - // Create variable to avoid multiple look-up - const bool prodTargetVisible = geomOptions->isVisible( "ProductionTarget" ); - const bool prodTargetSolid = geomOptions->isSolid ( "ProductionTarget" ); - const bool forceAuxEdgeVisible = geomOptions->forceAuxEdgeVisible( "ProductionTarget" ); - const bool placePV = geomOptions->placePV( "ProductionTarget" ); - const bool doSurfaceCheck = geomOptions->doSurfaceCheck( "ProductionTarget" ); - - for (std::map::const_iterator iPnt=anchoringPntsRgt.begin(); iPnt!=anchoringPntsRgt.end(); ++iPnt) { - double tmpAngle = iPnt->first * CLHEP::deg; - CLHEP::Hep3Vector normalax(cos(tmpAngle),sin(tmpAngle),0.0); - CLHEP::Hep3Vector tmpAnchPnt = _supWheel_trgtPS_rIn*normalax; - CLHEP::Hep3Vector tmpDirVec = tmpAnchPnt-iPnt->second; - CLHEP::Hep3Vector tmpMidPnt(0.5*(tmpAnchPnt.x()+iPnt->second.x()),0.5*(tmpAnchPnt.y()+iPnt->second.y()),0.5*(tmpAnchPnt.z()+iPnt->second.z())); - double tmpSpokeLength = tmpDirVec.mag(); - tmpDirVec = tmpDirVec.unit(); - CLHEP::Hep3Vector tmpRotAxis = tmpDirVec.cross(zax); - double rotAngle = tmpDirVec.angle(zax); - double deepAngleWheel = tmpDirVec.angle(normalax); - - CLHEP::Hep3Vector normalaxHub(sin(normalOverhangRgt),0,cos(normalOverhangRgt)); - normalaxHub.rotateZ(tmpAngle); - normalaxHub.transform(invTrgtRot); - double deepAngleHub = tmpDirVec.angle(normalaxHub); - - double fixOverlapWheel = spokeRad*tan(deepAngleWheel); - double fixOverlapHub = spokeRad*tan(deepAngleHub); - tmpSpokeLength-=(fixOverlapWheel+fixOverlapHub); - tmpMidPnt += 0.5*(fixOverlapHub-fixOverlapWheel)*tmpDirVec; - - TubsParams iSpokeParams( 0.0, spokeRad, 0.5*tmpSpokeLength); - std::stringstream iSpokeName; - iSpokeName<<"ProductionTargetSpokeRgt_"<startingSectionThickness().at(ithSection)/2.); + double currentHalfStartingSegment = tgt->startingSectionThickness().at(ithSection)/2.; + // + //set currentZ to be in the center of the starting section + _currentZ += tgt->startingSectionThickness().at(ithSection)/2.; + if (verbosityLevel > 3){ + G4cout << __PRETTY_FUNCTION__ << " ithSection , startingSegmentCenter at " << ithSection << " " <<_currentZ << G4endl; + G4cout << " and copy number for starting segment is now " << coreCopyNumber << G4endl; } + _startingSegmentCenter.setZ(_currentZ); - std::map const & anchoringPntsLft = tgt->anchoringPntsLft(); - - iSpk=0; - for (std::map::const_iterator iPnt=anchoringPntsLft.begin(); iPnt!=anchoringPntsLft.end(); ++iPnt) { - double tmpAngle = iPnt->first * CLHEP::deg; - CLHEP::Hep3Vector normalax(cos(tmpAngle),sin(tmpAngle),0.0); - CLHEP::Hep3Vector tmpAnchPnt = _supWheel_trgtPS_rIn*normalax; - CLHEP::Hep3Vector tmpDirVec = tmpAnchPnt-iPnt->second; - CLHEP::Hep3Vector tmpMidPnt(0.5*(tmpAnchPnt.x()+iPnt->second.x()),0.5*(tmpAnchPnt.y()+iPnt->second.y()),0.5*(tmpAnchPnt.z()+iPnt->second.z())); - double tmpSpokeLength = tmpDirVec.mag(); - tmpDirVec = tmpDirVec.unit(); - CLHEP::Hep3Vector tmpRotAxis = tmpDirVec.cross(zax); - double rotAngle = tmpDirVec.angle(zax); - double deepAngleWheel = tmpDirVec.angle(normalax); - - CLHEP::Hep3Vector normalaxHub(sin(normalOverhangLft),0,cos(normalOverhangLft)); - normalaxHub.rotateZ(tmpAngle); - normalaxHub.transform(invTrgtRot); - double deepAngleHub = tmpDirVec.angle(normalaxHub); - //CLHEP::HepRotation *tmpSpokeRot = reg.add(CLHEP::HepRotation(tmpRotAxis,rotAngle)); - - double fixOverlapWheel = spokeRad*tan(deepAngleWheel); - double fixOverlapHub = spokeRad*tan(deepAngleHub); - tmpSpokeLength-=(fixOverlapWheel+fixOverlapHub); - tmpMidPnt += 0.5*(fixOverlapHub-fixOverlapWheel)*tmpDirVec; - - TubsParams iSpokeParams( 0.0, spokeRad, 0.5*tmpSpokeLength); - std::stringstream iSpokeName; - iSpokeName<<"ProductionTargetSpokeLft_"<productionTargetRotation().inverse()*_startingSegmentCenter; + VolumeInfo startingSegment = nestTubs(startName, + startingSegmentParams, + prodTargetCoreMaterial, + &tgt->productionTargetRotation(), + currentStartingSegmentCenter, prodTargetMotherInfo, - 0, + ithSection, prodTargetVisible, - G4Colour::Gray(), + G4Colour::Yellow(), prodTargetSolid, forceAuxEdgeVisible, placePV, doSurfaceCheck ); - ++iSpk; - } - } - else if (_config.getInt("targetPS_version") == ProductionTargetMaker::hayman_v_2_0){ - int verbosityLevel = _config.getInt("PSHayman.verbosityLevel"); - verbosityLevel >0 && - cout << __func__ << " verbosityLevel on Hayman2.0 : " << verbosityLevel << endl; - G4ThreeVector _hallOriginInMu2e = parent.centerInMu2e(); - // Create variable to avoid multiple look-up - - G4GeometryOptions* geomOptions = art::ServiceHandle()->geomOptions(); - - bool prodTargetVisible = geomOptions->isVisible( "ProductionTarget" ); - bool prodTargetSolid = geomOptions->isSolid ( "ProductionTarget" ); - bool forceAuxEdgeVisible = geomOptions->forceAuxEdgeVisible( "ProductionTarget" ); - bool placePV = geomOptions->placePV( "ProductionTarget" ); - bool doSurfaceCheck = geomOptions->doSurfaceCheck( "ProductionTarget" ); - // + //build fins around starting segment + for (int ithFin = 0; ithFin < tgt->nHaymanFins(); ++ithFin) + { + const std::string name = "ProductionTargetFinStartingSection_" + std::to_string(ithSection+1) + "_Fin_" + std::to_string(ithFin); + if (verbosityLevel > 1) + {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngles.at(ithFin) << " fin copy number = " << finCopyNumber << G4endl;} - // begin all names with ProductionTarget so when we build sensitive detectors we can make them all - // sensitive at once with LVname.find("ProductionTarget") !=std::string::npos - // - // Build the production target. - GeomHandle tgt; - G4Material* prodTargetCoreMaterial = findMaterialOrThrow(tgt->targetCoreMaterial()); - G4Material* prodTargetFinMaterial = findMaterialOrThrow(tgt->targetFinMaterial()); - G4Material* prodTargetSupportRingMaterial = findMaterialOrThrow(tgt->supportRingMaterial()); - - TubsParams prodTargetMotherParams( 0. - ,tgt->productionTargetMotherOuterRadius() - ,tgt->productionTargetMotherHalfLength()); - - G4ThreeVector _loclCenter(0.0,0.0,0.0); - G4ThreeVector zeroTranslation(0.,0.,0.); - G4RotationMatrix* targetRotation = reg.add(G4RotationMatrix(tgt->productionTargetRotation().inverse())); - if (verbosityLevel > 2){G4cout << __PRETTY_FUNCTION__ << "target rotation = " << *targetRotation << G4endl;} - VolumeInfo prodTargetMotherInfo = nestTubs( "ProductionTargetMother", - prodTargetMotherParams, - parent.logical->GetMaterial(), - 0, - tgt->haymanProdTargetPosition() - parent.centerInMu2e(), - parent, - 0, - G4Colour::Blue(), - "PS" - ); + double finHalfHeightAboveTarget = (tgt->finOuterRadius() - tgt->rOut())/2.; + G4Box* finBox = new G4Box("finBox",tgt->haymanFinThickness()/2.,finHalfHeightAboveTarget,currentHalfStartingSegment); - if (verbosityLevel > 0){ - G4cout << "target position and hall origin = " - << tgt->haymanProdTargetPosition() << "\n" << - _hallOriginInMu2e << " " << parent.centerInMu2e() << G4endl; - } - if (verbosityLevel > 2){ - G4cout << __PRETTY_FUNCTION__ << "created prodTargetMotherInfo " - << tgt->productionTargetMotherOuterRadius() - << " " <productionTargetMotherHalfLength() << G4endl; - } - // - // construct each section + // + // this tubs is in the xy plane with an extent along z. So phi goes in the xy plane. + // need extra length for visualization tool and overlap avoidance + G4Tubs* tubCutout = new G4Tubs("finCutout",0.,tgt->rOut(),currentHalfStartingSegment + magicOffset,dSphi,dPphi); + ++finCopyNumber; - int numberOfSections = tgt->numberOfTargetSections(); + G4ThreeVector finTranslation = currentStartingSegmentCenter; + // + // G4 translates and then rotates. I want to offset the fin to its final location relative to the core and then let G4 perform rotations, + // both of the fin about the z axis and then the entire production target rotation angle + double distanceFromCenter = finHalfHeightAboveTarget + tgt->rOut(); + CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*std::cos(currentFinAngles.at(ithFin)) + ,distanceFromCenter*std::sin(currentFinAngles.at(ithFin)),0.); + // + // finTranslation was already put in the target rotated frame. I have to rotate the offset vector as well + offsetRelativeToCore *= tgt->productionTargetRotation().inverse(); + finTranslation = finTranslation + offsetRelativeToCore; + + CLHEP::Hep3Vector downshift = CLHEP::Hep3Vector(0.,-finHalfHeightAboveTarget - targetRadius*std::cos(angularSize), 0.); + placeSubtractionVolumeBoxTubs(name + ,finTranslation + ,prodTargetMotherInfo + ,finBox + ,tubCutout + ,rotFinsG4.at(ithFin) + ,prodTargetFinMaterial + ,nullptr + ,downshift + ,finCopyNumber + ,G4Colour::Blue() + ,"PS" + ,verbosityLevel>1); + // + // for each fin, build the little section at the top that joins it to the next fin. This will be a G4Box with a cutout. + // Here make the cutout a little tubs without the fancy angular business, since I don't have fins overlapping fins up at the top. + const std::string nameTop = "ProductionTargetFinTopStartingSection_" + std::to_string(ithSection+1) + "_Fin_" + std::to_string(ithFin); + + double gapRadius = tgt->thicknessOfGapPerSection().at(ithSection)/2.; + double finTopHalfHeight =( tgt->finOuterRadius() - tgt->heightOfRectangularGapPerSection().at(ithSection))/2.; + if (verbosityLevel > 2){G4cout << "finTopHalfHeight = " << finTopHalfHeight <haymanFinThickness()/2.,finTopHalfHeight,tgt->thicknessOfGapPerSection().at(ithSection)/2.); + G4Tubs* finTopCutout = new G4Tubs("finTopCutout",0.,gapRadius + magicOffset,tgt->haymanFinThickness()/2. + magicOffset,0.,2.*M_PI); + ++finTopCopyNumber; + + if (verbosityLevel > 2){G4cout << "current starting SegmentCenter, current HalfSegment, gap radius = " + << currentStartingSegmentCenter << " " << currentHalfStartingSegment + << " " << gapRadius << G4endl;} + G4ThreeVector finTopLocation = currentStartingSegmentCenter + + tgt->productionTargetRotation().inverse()*CLHEP::Hep3Vector(0.,0.,currentHalfStartingSegment + gapRadius); + // + // for debugging + CLHEP::Hep3Vector offsetTop(0.,0.,0.); + // + // rotate the offset vector + G4ThreeVector finTopTranslation = finTopLocation + offsetTop; + double distanceTopFromCenter = tgt->finOuterRadius() - finTopHalfHeight; + CLHEP::Hep3Vector offsetTopRelativeToCore(distanceTopFromCenter*std::cos(currentFinAngles.at(ithFin)) + ,distanceTopFromCenter*std::sin(currentFinAngles.at(ithFin)),0.); + // + // finTranslation was already put in the target rotated frame. I have to rotate the offset vector as well + offsetTopRelativeToCore *= tgt->productionTargetRotation().inverse(); + finTopTranslation = finTopTranslation + offsetTopRelativeToCore; + if (verbosityLevel > 2){G4cout << "finTopTranslation = " << finTopTranslation << G4endl;} + CLHEP::Hep3Vector downshiftTopBox = CLHEP::Hep3Vector(0.,-finTopHalfHeight, 0.); + // + // these are oriented with gap radius along z. Recall still in mother volume so this axis is OK + G4RotationMatrix* tubTopRotation = reg.add(G4RotationMatrix(CLHEP::HepRotation::IDENTITY)); + tubTopRotation->rotateY(M_PI/2.); - // - // first build the core + VolumeInfo finTopWithCutoutInfo(nameTop,finTopTranslation,prodTargetMotherInfo.centerInWorld); + finTopWithCutoutInfo.solid = new G4SubtractionSolid(nameTop,finTopBox,finTopCutout,tubTopRotation,downshiftTopBox); + finishNesting(finTopWithCutoutInfo + ,prodTargetFinMaterial + ,rotFinsG4.at(ithFin) + ,finTopTranslation + ,prodTargetMotherInfo.logical + ,finTopCopyNumber + ,G4Colour::White() + ,"ProductionTarget" + ,verbosityLevel>1); + // + // this check also uses finWithCutout. Indirectly it's a way of forcing an overlap check or it won't compile, since finWithCutout is never used otherwise. + // if (doSurfaceCheck){checkForOverlaps(finTopWithCutout,_config,0);} + + } // - // start at most negative end (technically could start at center but this is simpler since target is not symmetric about its center - // this variable gets incremented every time we add anything - double _currentZ = _loclCenter.z() - tgt->halfHaymanLength(); - if (verbosityLevel > 2){G4cout << __PRETTY_FUNCTION__ << " current Z starts at " << _currentZ << G4endl;} - // - // running z locations for placing elements - CLHEP::Hep3Vector _segmentCenter(0.,0.,0.); - CLHEP::Hep3Vector _startingSegmentCenter(0.,0.,0.); - double targetRadius = tgt->rOut(); + // for very first starting section, we need to build the support ring // - // keeping track of these for sensitive volume creation - int coreCopyNumber = 0; - int finCopyNumber = 0; - int finTopCopyNumber = 0; + // build support ring at this end. I need to build an annulus with cutouts for the fins. + int boxCopyNumber = 0; + if (ithSection==0){ + std::string name = "ProductionTargetNegativeEndRing"; + G4Tubs* supportRing = new G4Tubs(name, + tgt->supportRingInnerRadius(), + tgt->supportRingOuterRadius(), + tgt->supportRingLength()/2., + 0., + 2.*M_PI); + + if (verbosityLevel > 0){ + G4cout << __PRETTY_FUNCTION__ << ": \n " << name.c_str() << ":\n (rin, rout, half length) = (" + << tgt->supportRingInnerRadius() << ", " << tgt->supportRingOuterRadius() + << ", " << tgt->supportRingLength()/2. << ")" << G4endl; + } + // move ring to end of target and then back by cutout size (signs for negative side) + G4ThreeVector ringTranslation = currentStartingSegmentCenter + - tgt->productionTargetRotation().inverse()*CLHEP::Hep3Vector(0.,0.,tgt->supportRingLength()/2.); - double angularSize = asin((tgt->haymanFinThickness()/2.)/targetRadius); - double dSphi = M_PI/2. - angularSize; - double dPphi = 2.*angularSize; + if (verbosityLevel > 0){ + G4cout << " center = (0., 0., -" << (currentStartingSegmentCenter.mag() + tgt->supportRingLength()/2.) << ")" << G4endl; + G4cout << " center (wrt target mother) = " << ringTranslation << G4endl; + } + // the way we handle the cutout is a little different; the tubs is centered on zero, not offset from the z axis + // build one cutout box for each fin. This ignores the extra cutouts (easy enough to put in later) and the mounting mechanisms, + // which won't matter for muon yield, heat deposition, etc -- just looks like tungsten for those purposes - // - // need this for a) drawing subtraction volumes and b) rectangular cutouts out of circular support rings, etc. otherwise visualizer gets - // confused and how can G4 know which of two volumes some plane, for example, belongs to? - constexpr double magicOffset = 0.0001; + // + // this is all so I can add multiple cutouts to the ring and only do one placement. + // std::vector ringWithCutoutSolid; + std::vector ringWithCutoutSolid; - // - //other constants - double cutoutHeightAboveTarget = tgt->supportRingInnerRadius() + tgt->supportRingCutoutThickness()/2. - magicOffset; + std::string nameRing = "ProductionTargetNegativeRingCutout"; + + for (int ithFin = 0; ithFin < tgt->nHaymanFins(); ++ithFin){ + double currentFinAngle = tgt->finAngles().at(ithFin); + const std::string name = "ProductionTargetNegativeRingCutout_" + std::to_string(ithSection) + + "_Segment_" +"_Fin_" + std::to_string(ithFin); + ++boxCopyNumber; + + if (verbosityLevel > 6) + {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngle << " fin copy number = " << finCopyNumber << G4endl;} + if (verbosityLevel > 2){G4cout << " cutoutHeightAboveTarget = " << tgt->supportRingInnerRadius() << " " + << " " << tgt->supportRingCutoutThickness() << " " << cutoutHeightAboveTarget << G4endl;} + + // + // G4 translates and then rotates. I want to offset the fin to its final location relative to the core and then let G4 perform rotations, + // both of the fin about the z axis and then the entire production target rotation angle + double distanceFromCenter = cutoutHeightAboveTarget; + CLHEP::Hep3Vector boxShift(distanceFromCenter*std::cos(currentFinAngles.at(ithFin)),distanceFromCenter*std::sin(currentFinAngles.at(ithFin)),0.); + CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*std::cos(currentFinAngles.at(ithFin)),distanceFromCenter*std::sin(currentFinAngles.at(ithFin)), + tgt->supportRingCutoutLength()/2. + magicOffset); + + if (verbosityLevel > 3){ + G4cout << "production target rotation angle = " << tgt->productionTargetRotation().getTheta() << G4endl; + G4cout << "fin dump: " << currentFinAngles.at(ithFin) << " " << distanceFromCenter << " " << offsetRelativeToCore << G4endl; + G4cout << "ithfin; current fin angle; rotFin; finTranslation = " << ithFin << " " + << currentFinAngles.at(ithFin) << " " << *rotFinsG4.at(ithFin) << " " << ringTranslation << G4endl; + G4cout << "production target rotation = " << tgt->productionTargetRotation() << G4endl; + } + + if (ithFin == 0){ + ringWithCutoutSolid.push_back(new G4SubtractionSolid(name,supportRing,cutoutBox,rotBoxesG4.at(ithFin),offsetRelativeToCore)); + }else { + ringWithCutoutSolid.push_back(new G4SubtractionSolid(name,ringWithCutoutSolid.at(ithFin-1),cutoutBox,rotBoxesG4.at(ithFin),offsetRelativeToCore)); + } + } + + + if (verbosityLevel > 0){ + G4cout << " center (at end) = (0., 0., " << ringTranslation.mag() << ")" <1); - CLHEP::HepRotation* rotRing = reg.add(CLHEP::HepRotation(tgt->productionTargetRotation())); - // - // we're going to use these over and over; save time, memory, make cleaner and put in registry. save typing with currentFinAngles - // fins are constructed at 90^o relative to boxes and it's just cleaner to do it this way - std::vector rotFinsG4; - std::vector rotBoxesG4; - std::vector currentFinAngles; - for (int ithFin = 0; ithFin < tgt->nHaymanFins(); ++ithFin){ - rotFinsG4.emplace_back(reg.add(CLHEP::HepRotation(tgt->productionTargetRotation()))); - rotBoxesG4.emplace_back(reg.add(CLHEP::HepRotation::IDENTITY)); - // - // here I have a little confusing fix. The fin is built along the y-axis. But the rotation is given wrt the x axis. Hence I need to - // subtract off that 90^o . Then the - sign is CLHEP vs G4 - currentFinAngles.emplace_back(tgt->finAngles().at(ithFin)); - rotFinsG4.at(ithFin)->rotateZ(-(-M_PI/2. + currentFinAngles.at(ithFin))); - rotBoxesG4.at(ithFin)->rotateZ(-(-M_PI/2. + currentFinAngles.at(ithFin))); - if (verbosityLevel > 1){G4cout << "rotFin " << ithFin << " " << *rotFinsG4.at(ithFin) << G4endl;} } - G4Box* cutoutBox = reg.add(G4Box("cutoutBox",tgt->supportRingCutoutThickness()/2.+ magicOffset, - tgt->haymanFinThickness() + magicOffset,tgt->supportRingCutoutLength()/2.+ magicOffset)); + /* -------------------end of support ring, first starting section ---------------------*/ - // get real - for (int ithSection = 0; ithSection < numberOfSections; ++ithSection){ - int numberOfSegments = tgt->numberOfSegmentsPerSection().at(ithSection); - std::string startName = "ProductionTargetStartingCoreSection_" + std::to_string(ithSection+1); - // - // first place starting segment + // + // set current z to be at the end of this starting segment before beginning loop on sections + _currentZ += tgt->startingSectionThickness().at(ithSection)/2.; + if (verbosityLevel > 0){G4cout << __PRETTY_FUNCTION__ << " ithSection , startingSegment Z ends at " << ithSection << " " <<_currentZ << G4endl;} + double currentGap = tgt->thicknessOfGapPerSection().at(ithSection); + double currentHalfSegment = tgt->thicknessOfSegmentPerSection().at(ithSection)/2.; + + for (int ithSegment = 0; ithSegment < numberOfSegments; ++ithSegment){ + std::string name = "ProductionTargetCoreSection_" + std::to_string(ithSection+1) + + "_Segment_" + std::to_string(ithSegment); + if (verbosityLevel > 0) { + G4cout << __PRETTY_FUNCTION__ << "name = " << name << " and core copy number = " << coreCopyNumber << G4endl; + G4cout << __PRETTY_FUNCTION__ << "gap, and segment half = " << currentGap << " " << currentHalfSegment << G4endl; + } + TubsParams segmentParams(0.,targetRadius,currentHalfSegment); + _currentZ += currentHalfSegment + currentGap; + if (verbosityLevel > 0){ + G4cout << __PRETTY_FUNCTION__ << " ithSection , current Z for segment center is at " << ithSection << " " + << ithSegment << " " <<_currentZ << G4endl; + } + _segmentCenter.setZ(_currentZ); + G4ThreeVector currentSegmentCenter = tgt->productionTargetRotation().inverse()*_segmentCenter; + //CLHEP::Hep3Vector currentSegmentCenter = tgt->productionTargetRotation()*_segmentCenter; + // CLHEP::Hep3Vector currentSegmentCenter = _segmentCenter; + VolumeInfo coreSegment = nestTubs(name, + segmentParams, + prodTargetCoreMaterial, + &tgt->productionTargetRotation(), + currentSegmentCenter, + prodTargetMotherInfo, + coreCopyNumber, + prodTargetVisible, + G4Colour::Yellow(), + prodTargetSolid, + forceAuxEdgeVisible, + placePV, + doSurfaceCheck + ); + if (verbosityLevel > 2){G4cout << "segment center after volumeinfo = " << _segmentCenter << G4endl;} - TubsParams startingSegmentParams(0.,targetRadius,tgt->startingSectionThickness().at(ithSection)/2.); - double currentHalfStartingSegment = tgt->startingSectionThickness().at(ithSection)/2.; // - //set currentZ to be in the center of the starting section - _currentZ += tgt->startingSectionThickness().at(ithSection)/2.; - if (verbosityLevel > 3){ - G4cout << __PRETTY_FUNCTION__ << " ithSection , startingSegmentCenter at " << ithSection << " " <<_currentZ << G4endl; - G4cout << " and copy number for starting segment is now " << coreCopyNumber << G4endl; - } - _startingSegmentCenter.setZ(_currentZ); + //now add fins surrounding this core segment. Build them as simple rectangles with a G4 subtraction volume for the core. - CLHEP::Hep3Vector currentStartingSegmentCenter = tgt->productionTargetRotation().inverse()*_startingSegmentCenter; - VolumeInfo startingSegment = nestTubs(startName, - startingSegmentParams, - prodTargetCoreMaterial, - &tgt->productionTargetRotation(), - currentStartingSegmentCenter, - prodTargetMotherInfo, - ithSection, - prodTargetVisible, - G4Colour::Yellow(), - prodTargetSolid, - forceAuxEdgeVisible, - placePV, - doSurfaceCheck - ); - //build fins around starting segment for (int ithFin = 0; ithFin < tgt->nHaymanFins(); ++ithFin) + // for (int ithFin = 0; ithFin < 1; ++ithFin) { - const std::string name = "ProductionTargetFinStartingSection_" + std::to_string(ithSection+1) + "_Fin_" + std::to_string(ithFin); + const std::string name = "ProductionTargetFinSection_" + std::to_string(ithSection+1) + "_Segment_" + std::to_string(ithSegment) + "_Fin_" + std::to_string(ithFin); if (verbosityLevel > 1) {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngles.at(ithFin) << " fin copy number = " << finCopyNumber << G4endl;} - + // divide by 2 to make box half-height since I gave it the "radius" to start with; arbitrary convention. double finHalfHeightAboveTarget = (tgt->finOuterRadius() - tgt->rOut())/2.; - G4Box* finBox = new G4Box("finBox",tgt->haymanFinThickness()/2.,finHalfHeightAboveTarget,currentHalfStartingSegment); - + if (verbosityLevel > 2){G4cout << "finHeightAboveTarget = " << tgt->finOuterRadius() << " " + << tgt->rOut() << " " << finHalfHeightAboveTarget << G4endl;} + G4Box* finBox = new G4Box("finBox",tgt->haymanFinThickness()/2.,finHalfHeightAboveTarget,currentHalfSegment); // // this tubs is in the xy plane with an extent along z. So phi goes in the xy plane. - // need extra length for visualization tool and overlap avoidance - G4Tubs* tubCutout = new G4Tubs("finCutout",0.,tgt->rOut(),currentHalfStartingSegment + magicOffset,dSphi,dPphi); + G4Tubs* tubCutout = new G4Tubs("finCutout",0.,tgt->rOut(),currentHalfSegment + magicOffset,dSphi,dPphi);// need extra length for visualization tool + // + // I now have two G4Solids, a box and a cutout. Combine them to make a logical volume. the box is + // oriented so that it is "z" thick and "radius" high. the tubs is made to be the same, so no rotation + // and also no translation needed + ++finCopyNumber; - G4ThreeVector finTranslation = currentStartingSegmentCenter; // - // G4 translates and then rotates. I want to offset the fin to its final location relative to the core and then let G4 perform rotations, - // both of the fin about the z axis and then the entire production target rotation angle + //for debugging + CLHEP::Hep3Vector offset(0.,0.,0.); + G4ThreeVector finTranslation = currentSegmentCenter; + double distanceFromCenter = finHalfHeightAboveTarget + tgt->rOut(); - CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*cos(currentFinAngles.at(ithFin)) - ,distanceFromCenter*sin(currentFinAngles.at(ithFin)),0.); + CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*std::cos(currentFinAngles.at(ithFin)),distanceFromCenter*std::sin(currentFinAngles.at(ithFin)),0.); // // finTranslation was already put in the target rotated frame. I have to rotate the offset vector as well offsetRelativeToCore *= tgt->productionTargetRotation().inverse(); finTranslation = finTranslation + offsetRelativeToCore; - CLHEP::Hep3Vector downshift = CLHEP::Hep3Vector(0.,-finHalfHeightAboveTarget - targetRadius*cos(angularSize), 0.); + + CLHEP::Hep3Vector downshift = CLHEP::Hep3Vector(0.,-finHalfHeightAboveTarget - targetRadius*std::cos(angularSize), 0.); placeSubtractionVolumeBoxTubs(name ,finTranslation ,prodTargetMotherInfo @@ -609,646 +880,1427 @@ namespace mu2e { ,G4Colour::Blue() ,"PS" ,verbosityLevel>1); + // + // for each fin, build the little section at the top that joins it to the next fin. This will be a G4Box with a cutout. + // Here make the cutout a little tubs without the fancy angular business, since I don't have fins overlapping fins up at the top. + const std::string nameTop = "ProductionTargetFinTopSection_" + std::to_string(ithSection+1) + + "_Segment_" + std::to_string(ithSegment) + "_Fin_" + std::to_string(ithFin); + + //Much easier to visualize but the code is equally yucky + double gapRadius = tgt->thicknessOfGapPerSection().at(ithSection)/2.; + double finTopHalfHeight =( tgt->finOuterRadius() - tgt->heightOfRectangularGapPerSection().at(ithSection))/2.; + if (verbosityLevel > 2){G4cout << "finTopHalfHeight = " << finTopHalfHeight <haymanFinThickness()/2.,finTopHalfHeight,tgt->thicknessOfGapPerSection().at(ithSection)/2.); + G4Tubs* finTopCutout = new G4Tubs("finTopCutout",0.,gapRadius + magicOffset,tgt->haymanFinThickness()/2. + magicOffset,0.,2.*M_PI); + ++finTopCopyNumber; + + if (verbosityLevel > 2){G4cout << "current starting SegmentCenter, current HalfSegment, gap radius = " + << currentStartingSegmentCenter << " " << currentHalfStartingSegment << " " + << gapRadius << G4endl;} + G4ThreeVector finTopBoxLocation = currentSegmentCenter + + tgt->productionTargetRotation().inverse()*CLHEP::Hep3Vector(0.,0.,currentHalfSegment + gapRadius); + // + // for debugging + CLHEP::Hep3Vector offsetTop(0.,0.,0.); + // + // rotate the offset vector + G4ThreeVector finTopTranslation = finTopBoxLocation + offsetTop; + // + // G4 translates and then rotates. I want to offset the fin to its final location relative to the core and then let G4 perform rotations, + // both of the fin about the z axis and then the entire production target rotation angle + double distanceTopFromCenter = tgt->finOuterRadius() - finTopHalfHeight; - // - // for each fin, build the little section at the top that joins it to the next fin. This will be a G4Box with a cutout. - // Here make the cutout a little tubs without the fancy angular business, since I don't have fins overlapping fins up at the top. - const std::string nameTop = "ProductionTargetFinTopStartingSection_" + std::to_string(ithSection+1) + "_Fin_" + std::to_string(ithFin); - - double gapRadius = tgt->thicknessOfGapPerSection().at(ithSection)/2.; - double finTopHalfHeight =( tgt->finOuterRadius() - tgt->heightOfRectangularGapPerSection().at(ithSection))/2.; - if (verbosityLevel > 2){G4cout << "finTopHalfHeight = " << finTopHalfHeight <haymanFinThickness()/2.,finTopHalfHeight,tgt->thicknessOfGapPerSection().at(ithSection)/2.); - G4Tubs* finTopCutout = new G4Tubs("finTopCutout",0.,gapRadius + magicOffset,tgt->haymanFinThickness()/2. + magicOffset,0.,2.*M_PI); - ++finTopCopyNumber; - - if (verbosityLevel > 2){G4cout << "current starting SegmentCenter, current HalfSegment, gap radius = " - << currentStartingSegmentCenter << " " << currentHalfStartingSegment - << " " << gapRadius << G4endl;} - G4ThreeVector finTopLocation = currentStartingSegmentCenter - + tgt->productionTargetRotation().inverse()*CLHEP::Hep3Vector(0.,0.,currentHalfStartingSegment + gapRadius); - // - // for debugging - CLHEP::Hep3Vector offsetTop(0.,0.,0.); - // - // rotate the offset vector - G4ThreeVector finTopTranslation = finTopLocation + offsetTop; - double distanceTopFromCenter = tgt->finOuterRadius() - finTopHalfHeight; - CLHEP::Hep3Vector offsetTopRelativeToCore(distanceTopFromCenter*cos(currentFinAngles.at(ithFin)) - ,distanceTopFromCenter*sin(currentFinAngles.at(ithFin)),0.); - // - // finTranslation was already put in the target rotated frame. I have to rotate the offset vector as well - offsetTopRelativeToCore *= tgt->productionTargetRotation().inverse(); - finTopTranslation = finTopTranslation + offsetTopRelativeToCore; - if (verbosityLevel > 2){G4cout << "finTopTranslation = " << finTopTranslation << G4endl;} - CLHEP::Hep3Vector downshiftTopBox = CLHEP::Hep3Vector(0.,-finTopHalfHeight, 0.); - // - // these are oriented with gap radius along z. Recall still in mother volume so this axis is OK - G4RotationMatrix* tubTopRotation = reg.add(G4RotationMatrix(CLHEP::HepRotation::IDENTITY)); - - tubTopRotation->rotateY(M_PI/2.); - - VolumeInfo finTopWithCutoutInfo(nameTop,finTopTranslation,prodTargetMotherInfo.centerInWorld); - finTopWithCutoutInfo.solid = new G4SubtractionSolid(nameTop,finTopBox,finTopCutout,tubTopRotation,downshiftTopBox); - finishNesting(finTopWithCutoutInfo - ,prodTargetFinMaterial - ,rotFinsG4.at(ithFin) - ,finTopTranslation - ,prodTargetMotherInfo.logical - ,finTopCopyNumber - ,G4Colour::White() - ,"ProductionTarget" - ,verbosityLevel>1); - // - // this check also uses finWithCutout. Indirectly it's a way of forcing an overlap check or it won't compile, since finWithCutout is never used otherwise. - // if (doSurfaceCheck){checkForOverlaps(finTopWithCutout,_config,0);} + CLHEP::Hep3Vector offsetTopRelativeToCore(distanceTopFromCenter*std::cos(currentFinAngles.at(ithFin)),distanceTopFromCenter*std::sin(currentFinAngles.at(ithFin)),0.); + // + // finTranslation was already put in the target rotated frame. I have to rotate the offset vector as well + offsetTopRelativeToCore *= tgt->productionTargetRotation().inverse(); + finTopTranslation = finTopTranslation + offsetTopRelativeToCore; + if (verbosityLevel > 2){G4cout << "finTopTranslation = " << finTopTranslation << G4endl;} + CLHEP::Hep3Vector downshiftTopBox = CLHEP::Hep3Vector(0.,-finTopHalfHeight, 0.); - } + // + // these are oriented with gap radius along z. Recall still in mother volume so this axis is OK + G4RotationMatrix* tubTopRotation = reg.add(G4RotationMatrix(CLHEP::HepRotation::IDENTITY)); + tubTopRotation->rotateY(M_PI/2.); + + VolumeInfo finTopWithCutoutInfo(nameTop,finTopTranslation,prodTargetMotherInfo.centerInWorld); + finTopWithCutoutInfo.solid = new G4SubtractionSolid(nameTop,finTopBox,finTopCutout,tubTopRotation,downshiftTopBox); + finishNesting(finTopWithCutoutInfo + ,prodTargetFinMaterial + ,rotFinsG4.at(ithFin) + ,finTopTranslation + ,prodTargetMotherInfo.logical + ,finTopCopyNumber + ,G4Colour::White() + ,"ProductionTarget" + ,verbosityLevel>1); + } + _currentZ += currentHalfSegment; // - // for very first starting section, we need to build the support ring + // if this is the last segment, add another gap before next section starts! + if (ithSegment == numberOfSegments-1) { + if (verbosityLevel > 0){G4cout << " z before = " << _currentZ << G4endl;} + _currentZ += currentGap; + if (verbosityLevel > 0){G4cout << " z after = " << _currentZ << G4endl;} + } + if (verbosityLevel > 0){G4cout << __PRETTY_FUNCTION__ << " ending at z=" << _currentZ << G4endl;} + } + ++coreCopyNumber; + // + // if this is the last section, add on one more beginning block. decided not to extend vectors and have zero segments, just ugly + if (ithSection == (numberOfSections -1) ){ + startName = "ProductionTargetStartingCoreSection_" + std::to_string(ithSection+1+1); + TubsParams startingSegmentParamsEnd(0.,targetRadius,tgt->startingSectionThickness().at(ithSection)/2.); + // - // build support ring at this end. I need to build an annulus with cutouts for the fins. - int boxCopyNumber = 0; - if (ithSection==0){ - std::string name = "ProductionTargetNegativeEndRing"; - G4Tubs* supportRing = new G4Tubs(name, - tgt->supportRingInnerRadius(), - tgt->supportRingOuterRadius(), - tgt->supportRingLength()/2., - 0., - 2.*M_PI); - - if (verbosityLevel > 0){ - G4cout << __PRETTY_FUNCTION__ << ": \n " << name.c_str() << ":\n (rin, rout, half length) = (" - << tgt->supportRingInnerRadius() << ", " << tgt->supportRingOuterRadius() - << ", " << tgt->supportRingLength()/2. << ")" << G4endl; - } - // move ring to end of target and then back by cutout size (signs for negative side) - G4ThreeVector ringTranslation = currentStartingSegmentCenter - - tgt->productionTargetRotation().inverse()*CLHEP::Hep3Vector(0.,0.,tgt->supportRingLength()/2.); + //set currentZ to be in the center of the starting section + _currentZ += tgt->startingSectionThickness().at(ithSection)/2.; + if (verbosityLevel > 0){G4cout << __PRETTY_FUNCTION__ << " ithSection+1 , startingSegmentCenter at " << ithSection+1 << " " <<_currentZ << G4endl;} + _startingSegmentCenter.setZ(_currentZ); + CLHEP::Hep3Vector currentStartingSegmentCenter = tgt->productionTargetRotation().inverse()*_startingSegmentCenter; + double currentHalfStartingSegment = tgt->startingSectionThickness().at(ithSection)/2.; + startingSegment = nestTubs(startName, + startingSegmentParamsEnd, + prodTargetCoreMaterial, + &tgt->productionTargetRotation(), + currentStartingSegmentCenter, + prodTargetMotherInfo, + ithSection+1, + prodTargetVisible, + G4Colour::Yellow(), + prodTargetSolid, + forceAuxEdgeVisible, + placePV, + doSurfaceCheck + ); - if (verbosityLevel > 0){ - G4cout << " center = (0., 0., -" << (currentStartingSegmentCenter.mag() + tgt->supportRingLength()/2.) << ")" << G4endl; - G4cout << " center (wrt target mother) = " << ringTranslation << G4endl; - } - // the way we handle the cutout is a little different; the tubs is centered on zero, not offset from the z axis - // build one cutout box for each fin. This ignores the extra cutouts (easy enough to put in later) and the mounting mechanisms, - // which won't matter for muon yield, heat deposition, etc -- just looks like tungsten for those purposes + //build fins around starting segment + for (int ithFin = 0; ithFin < tgt->nHaymanFins(); ++ithFin){ - // - // this is all so I can add multiple cutouts to the ring and only do one placement. - // std::vector ringWithCutoutSolid; - std::vector ringWithCutoutSolid; + double currentFinAngle = tgt->finAngles().at(ithFin); + const std::string name = "ProductionTargetFinStartingSection_" + std::to_string(ithSection+1+1) + "_Fin_" + std::to_string(ithFin); - std::string nameRing = "ProductionTargetNegativeRingCutout"; + if (verbosityLevel > 1) + {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngle << " fin copy number = " << finCopyNumber << G4endl;} - for (int ithFin = 0; ithFin < tgt->nHaymanFins(); ++ithFin){ - double currentFinAngle = tgt->finAngles().at(ithFin); - const std::string name = "ProductionTargetNegativeRingCutout_" + std::to_string(ithSection) - + "_Segment_" +"_Fin_" + std::to_string(ithFin); - ++boxCopyNumber; + double finHalfHeightAboveTarget = (tgt->finOuterRadius() - tgt->rOut())/2.; + G4Box* finBox = new G4Box("finBox",tgt->haymanFinThickness()/2.,finHalfHeightAboveTarget,currentHalfStartingSegment); + // + // this tubs is in the xy plane with an extent along z. So phi goes in the xy plane. + // need extra length for visualization tool and overlap avoidance + G4Tubs* tubCutout = new G4Tubs("finCutout",0.,tgt->rOut(),currentHalfStartingSegment + magicOffset,dSphi,dPphi); - if (verbosityLevel > 6) - {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngle << " fin copy number = " << finCopyNumber << G4endl;} - if (verbosityLevel > 2){G4cout << " cutoutHeightAboveTarget = " << tgt->supportRingInnerRadius() << " " - << " " << tgt->supportRingCutoutThickness() << " " << cutoutHeightAboveTarget << G4endl;} + CLHEP::HepRotation* rotFin = reg.add(CLHEP::HepRotation(tgt->productionTargetRotation().inverse())); + CLHEP::HepRotation* rotFinG4 = reg.add(CLHEP::HepRotation(tgt->productionTargetRotation())); + // + // ok here I have a little confusing fix. The fin is built along the y-axis. But the rotation is given wrt the x axis. Hence I need to + // subtract off that 90^o + double rotAdj = -M_PI/2. + currentFinAngle; + rotFin->rotateZ(rotAdj); + // + // g4 version. + rotFinG4->rotateZ(-rotAdj); + ++finCopyNumber; + // + // for debugging + CLHEP::Hep3Vector offset(0.,0.,0.); // - // G4 translates and then rotates. I want to offset the fin to its final location relative to the core and then let G4 perform rotations, - // both of the fin about the z axis and then the entire production target rotation angle - double distanceFromCenter = cutoutHeightAboveTarget; - CLHEP::Hep3Vector boxShift(distanceFromCenter*cos(currentFinAngles.at(ithFin)),distanceFromCenter*sin(currentFinAngles.at(ithFin)),0.); - CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*cos(currentFinAngles.at(ithFin)),distanceFromCenter*sin(currentFinAngles.at(ithFin)), - tgt->supportRingCutoutLength()/2. + magicOffset); - - if (verbosityLevel > 3){ - G4cout << "production target rotation angle = " << tgt->productionTargetRotation().getTheta() << G4endl; - G4cout << "fin dump: " << currentFinAngles.at(ithFin) << " " << distanceFromCenter << " " << offsetRelativeToCore << G4endl; - G4cout << "ithfin; current fin angle; rotFin; finTranslation = " << ithFin << " " - << currentFinAngles.at(ithFin) << " " << *rotFinsG4.at(ithFin) << " " << ringTranslation << G4endl; - G4cout << "production target rotation = " << tgt->productionTargetRotation() << G4endl; - } - - if (ithFin == 0){ - ringWithCutoutSolid.push_back(new G4SubtractionSolid(name,supportRing,cutoutBox,rotBoxesG4.at(ithFin),offsetRelativeToCore)); - }else { - ringWithCutoutSolid.push_back(new G4SubtractionSolid(name,ringWithCutoutSolid.at(ithFin-1),cutoutBox,rotBoxesG4.at(ithFin),offsetRelativeToCore)); - } - } - - - if (verbosityLevel > 0){ - G4cout << " center (at end) = (0., 0., " << ringTranslation.mag() << ")" <1); + } + // + // set current z to be at the end of this starting segment as a final check + _currentZ += tgt->startingSectionThickness().at(ithSection)/2.; + if (verbosityLevel > 0){ + G4cout << __PRETTY_FUNCTION__ << " ithSection+1 , startingSegment Z ends at " << ithSection+1 << " " <<_currentZ << G4endl; + G4cout << " with copy number for last segment = " << ithSection+1 << G4endl; + } + /***********and now the final end ring ***********/ + std::string name = "ProductionTargetPositiveEndRing"; + G4Tubs* supportRing = new G4Tubs(name, + tgt->supportRingInnerRadius(), + tgt->supportRingOuterRadius(), + tgt->supportRingLength()/2., + 0., + 2.*M_PI); + if (verbosityLevel > 0){ + G4cout << __PRETTY_FUNCTION__ << ": \n " << name.c_str() << ":\n (rin, rout, half length) = (" + << tgt->supportRingInnerRadius() << ", " << tgt->supportRingOuterRadius() + << ", " << tgt->supportRingLength()/2. << ")" << G4endl; + } + // move ring to end of target and then back by cutout size (signs for positive side) + G4ThreeVector ringTranslation = currentStartingSegmentCenter + + tgt->productionTargetRotation().inverse()*CLHEP::Hep3Vector(0.,0.,tgt->supportRingLength()/2.); + if (verbosityLevel > 0){ + G4cout << " center = (0., 0., " << (currentStartingSegmentCenter.mag() + tgt->supportRingLength()/2.) << ")" << G4endl; + G4cout << " center (wrt target mother) = " << ringTranslation << G4endl; } + std::vector ringWithCutoutSolid; + std::string nameRing = "ProductionTargetPositiveRingCutout"; + for (int ithFin = 0; ithFin < tgt->nHaymanFins(); ++ithFin){ + double currentFinAngle = tgt->finAngles().at(ithFin); + const std::string name = "ProductionTargetPositiveRingCutout_" + std::to_string(ithSection) + + "_Segment_" +"_Fin_" + std::to_string(ithFin); + ++boxCopyNumber; - /* -------------------end of support ring, first starting section ---------------------*/ - // - // set current z to be at the end of this starting segment before beginning loop on sections - _currentZ += tgt->startingSectionThickness().at(ithSection)/2.; - if (verbosityLevel > 0){G4cout << __PRETTY_FUNCTION__ << " ithSection , startingSegment Z ends at " << ithSection << " " <<_currentZ << G4endl;} - double currentGap = tgt->thicknessOfGapPerSection().at(ithSection); - double currentHalfSegment = tgt->thicknessOfSegmentPerSection().at(ithSection)/2.; - - for (int ithSegment = 0; ithSegment < numberOfSegments; ++ithSegment){ - std::string name = "ProductionTargetCoreSection_" + std::to_string(ithSection+1) - + "_Segment_" + std::to_string(ithSegment); - if (verbosityLevel > 0) { - G4cout << __PRETTY_FUNCTION__ << "name = " << name << " and core copy number = " << coreCopyNumber << G4endl; - G4cout << __PRETTY_FUNCTION__ << "gap, and segment half = " << currentGap << " " << currentHalfSegment << G4endl; - } - TubsParams segmentParams(0.,targetRadius,currentHalfSegment); - _currentZ += currentHalfSegment + currentGap; - if (verbosityLevel > 0){ - G4cout << __PRETTY_FUNCTION__ << " ithSection , current Z for segment center is at " << ithSection << " " - << ithSegment << " " <<_currentZ << G4endl; + // + // G4 translates and then rotates. I want to offset the fin to its final location relative to the core and then let G4 perform rotations, + // both of the fin about the z axis and then the entire production target rotation angle + double distanceFromCenter = cutoutHeightAboveTarget; + + CLHEP::Hep3Vector boxShift(distanceFromCenter*std::cos(currentFinAngle),distanceFromCenter*std::sin(currentFinAngle),0.); + CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*std::cos(currentFinAngle),distanceFromCenter*std::sin(currentFinAngle), + -tgt->supportRingCutoutLength()/2. - magicOffset);// note flipped sign from other end of target + + CLHEP::Hep3Vector downshift = CLHEP::Hep3Vector(0.,0.,0.);//for debugging + offsetRelativeToCore = offsetRelativeToCore + downshift; + if (ithFin == 0){ + ringWithCutoutSolid.push_back(new G4SubtractionSolid(name,supportRing,cutoutBox,rotBoxesG4.at(ithFin),offsetRelativeToCore)); + }else { + ringWithCutoutSolid.push_back(new G4SubtractionSolid(name,ringWithCutoutSolid.at(ithFin-1),cutoutBox,rotBoxesG4.at(ithFin),offsetRelativeToCore)); } - _segmentCenter.setZ(_currentZ); - G4ThreeVector currentSegmentCenter = tgt->productionTargetRotation().inverse()*_segmentCenter; - //CLHEP::Hep3Vector currentSegmentCenter = tgt->productionTargetRotation()*_segmentCenter; - // CLHEP::Hep3Vector currentSegmentCenter = _segmentCenter; - VolumeInfo coreSegment = nestTubs(name, - segmentParams, - prodTargetCoreMaterial, - &tgt->productionTargetRotation(), - currentSegmentCenter, - prodTargetMotherInfo, - coreCopyNumber, - prodTargetVisible, - G4Colour::Yellow(), - prodTargetSolid, - forceAuxEdgeVisible, - placePV, - doSurfaceCheck - ); - if (verbosityLevel > 2){G4cout << "segment center after volumeinfo = " << _segmentCenter << G4endl;} + } - // - //now add fins surrounding this core segment. Build them as simple rectangles with a G4 subtraction volume for the core. - - for (int ithFin = 0; ithFin < tgt->nHaymanFins(); ++ithFin) - // for (int ithFin = 0; ithFin < 1; ++ithFin) - { - const std::string name = "ProductionTargetFinSection_" + std::to_string(ithSection+1) + "_Segment_" + std::to_string(ithSegment) + "_Fin_" + std::to_string(ithFin); - if (verbosityLevel > 1) - {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngles.at(ithFin) << " fin copy number = " << finCopyNumber << G4endl;} - // divide by 2 to make box half-height since I gave it the "radius" to start with; arbitrary convention. - double finHalfHeightAboveTarget = (tgt->finOuterRadius() - tgt->rOut())/2.; - if (verbosityLevel > 2){G4cout << "finHeightAboveTarget = " << tgt->finOuterRadius() << " " - << tgt->rOut() << " " << finHalfHeightAboveTarget << G4endl;} - G4Box* finBox = new G4Box("finBox",tgt->haymanFinThickness()/2.,finHalfHeightAboveTarget,currentHalfSegment); - // - // this tubs is in the xy plane with an extent along z. So phi goes in the xy plane. - G4Tubs* tubCutout = new G4Tubs("finCutout",0.,tgt->rOut(),currentHalfSegment + magicOffset,dSphi,dPphi);// need extra length for visualization tool - // - // I now have two G4Solids, a box and a cutout. Combine them to make a logical volume. the box is - // oriented so that it is "z" thick and "radius" high. the tubs is made to be the same, so no rotation - // and also no translation needed - - ++finCopyNumber; - - // - //for debugging - CLHEP::Hep3Vector offset(0.,0.,0.); - G4ThreeVector finTranslation = currentSegmentCenter; - - double distanceFromCenter = finHalfHeightAboveTarget + tgt->rOut(); - CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*cos(currentFinAngles.at(ithFin)),distanceFromCenter*sin(currentFinAngles.at(ithFin)),0.); - // - // finTranslation was already put in the target rotated frame. I have to rotate the offset vector as well - offsetRelativeToCore *= tgt->productionTargetRotation().inverse(); - finTranslation = finTranslation + offsetRelativeToCore; - - - CLHEP::Hep3Vector downshift = CLHEP::Hep3Vector(0.,-finHalfHeightAboveTarget - targetRadius*cos(angularSize), 0.); - placeSubtractionVolumeBoxTubs(name - ,finTranslation - ,prodTargetMotherInfo - ,finBox - ,tubCutout - ,rotFinsG4.at(ithFin) - ,prodTargetFinMaterial - ,nullptr - ,downshift - ,finCopyNumber - ,G4Colour::Blue() - ,"PS" - ,verbosityLevel>1); - // - // for each fin, build the little section at the top that joins it to the next fin. This will be a G4Box with a cutout. - // Here make the cutout a little tubs without the fancy angular business, since I don't have fins overlapping fins up at the top. - const std::string nameTop = "ProductionTargetFinTopSection_" + std::to_string(ithSection+1) - + "_Segment_" + std::to_string(ithSegment) + "_Fin_" + std::to_string(ithFin); - - //Much easier to visualize but the code is equally yucky - double gapRadius = tgt->thicknessOfGapPerSection().at(ithSection)/2.; - double finTopHalfHeight =( tgt->finOuterRadius() - tgt->heightOfRectangularGapPerSection().at(ithSection))/2.; - if (verbosityLevel > 2){G4cout << "finTopHalfHeight = " << finTopHalfHeight <haymanFinThickness()/2.,finTopHalfHeight,tgt->thicknessOfGapPerSection().at(ithSection)/2.); - G4Tubs* finTopCutout = new G4Tubs("finTopCutout",0.,gapRadius + magicOffset,tgt->haymanFinThickness()/2. + magicOffset,0.,2.*M_PI); - ++finTopCopyNumber; - - if (verbosityLevel > 2){G4cout << "current starting SegmentCenter, current HalfSegment, gap radius = " - << currentStartingSegmentCenter << " " << currentHalfStartingSegment << " " - << gapRadius << G4endl;} - G4ThreeVector finTopBoxLocation = currentSegmentCenter - + tgt->productionTargetRotation().inverse()*CLHEP::Hep3Vector(0.,0.,currentHalfSegment + gapRadius); - // - // for debugging - CLHEP::Hep3Vector offsetTop(0.,0.,0.); - // - // rotate the offset vector - G4ThreeVector finTopTranslation = finTopBoxLocation + offsetTop; - // - // G4 translates and then rotates. I want to offset the fin to its final location relative to the core and then let G4 perform rotations, - // both of the fin about the z axis and then the entire production target rotation angle - double distanceTopFromCenter = tgt->finOuterRadius() - finTopHalfHeight; - CLHEP::Hep3Vector offsetTopRelativeToCore(distanceTopFromCenter*cos(currentFinAngles.at(ithFin)),distanceTopFromCenter*sin(currentFinAngles.at(ithFin)),0.); - // - // finTranslation was already put in the target rotated frame. I have to rotate the offset vector as well - offsetTopRelativeToCore *= tgt->productionTargetRotation().inverse(); - finTopTranslation = finTopTranslation + offsetTopRelativeToCore; - if (verbosityLevel > 2){G4cout << "finTopTranslation = " << finTopTranslation << G4endl;} - CLHEP::Hep3Vector downshiftTopBox = CLHEP::Hep3Vector(0.,-finTopHalfHeight, 0.); + if (verbosityLevel > 0){ + G4cout << " center (at end) = (0., 0., " << ringTranslation.mag() << ")" <1); - - } - _currentZ += currentHalfSegment; - // - // if this is the last segment, add another gap before next section starts! - if (ithSegment == numberOfSegments-1) { - if (verbosityLevel > 0){G4cout << " z before = " << _currentZ << G4endl;} - _currentZ += currentGap; - if (verbosityLevel > 0){G4cout << " z after = " << _currentZ << G4endl;} + G4cout << " rotation = " << rotRing << G4endl; + } + VolumeInfo ringWithCutoutPositive(name,ringTranslation,prodTargetMotherInfo.centerInWorld); + ringWithCutoutPositive.solid = ringWithCutoutSolid.at(tgt->nHaymanFins() - 1); // what does this = mean? + finishNesting(ringWithCutoutPositive + ,prodTargetSupportRingMaterial + ,rotRing + ,ringTranslation + ,prodTargetMotherInfo.logical + ,finCopyNumber + ,G4Colour::White() + ,"ProductionTarget" + ,verbosityLevel>1); + + } //end if(ithSection == (numberOfSections -1) ) + } //end for(int ithSection...) + + //Add support structures for the production target + if(tgt->supportsBuild()) { + G4Material* suppWheelMaterial = findMaterialOrThrow(tgt->supportWheelMaterial()); + G4ThreeVector localWheelCenter(0.0,0.0,0.0); //no offset + double suppWheelParams[] = {tgt->supportWheelRIn(), tgt->supportWheelROut(), tgt->supportWheelHL()}; + //create the volume info for the support wheel+rods + VolumeInfo suppWheelInfo( "ProductionTargetSupportWheel", localWheelCenter, prodTargetMotherInfo.centerInMu2e()); + suppWheelInfo.solid = new G4Tubs("ProductionTargetSupportWheel_wheel", suppWheelParams[0], suppWheelParams[1], + suppWheelParams[2], 0., CLHEP::twopi); + // suppWheelParams, + // suppWheelMaterial, + // 0, + // localWheelCenter, + // prodTargetMotherInfo, + // 0, + // G4Colour::Gray(), + // "PS" + // ); + + // add spokes // + + //spoke info + const int nspokesperside = tgt->nSpokesPerSide(); + G4Material* spokeMaterial = findMaterialOrThrow(tgt->spokeMaterial()); + //target info + double rTarget = tgt->supportRingOuterRadius(); //radius of the support ring to attach to + double zTarget = tgt->halfHaymanLength(); //where along the target to attach + double smallGap = 0.001; //for adding small offsets to avoid overlaps due to precision + //initialize parameter vectors + //features on wheel + const vector supportWheelFeatureAngles = tgt->supportWheelFeatureAngles(); + const vector supportWheelFeatureArcs = tgt->supportWheelFeatureArcs (); + const vector supportWheelFeatureRIns = tgt->supportWheelFeatureRIns (); + //support rods in wheel + const vector supportWheelRodHL = tgt->supportWheelRodHL (); + const vector supportWheelRodOffset = tgt->supportWheelRodOffset (); + const vector supportWheelRodRadius = tgt->supportWheelRodRadius (); + const vector supportWheelRodRadialOffset = tgt->supportWheelRodRadialOffset(); + const vector supportWheelRodWireOffsetD = tgt->supportWheelRodWireOffsetD (); + const vector supportWheelRodWireOffsetU = tgt->supportWheelRodWireOffsetU (); + const vector supportWheelRodAngles = tgt->supportWheelRodAngles (); + //spoke (support wire) angles + const vector spokeTargetAnglesD = tgt->spokeTargetAnglesD(); + const vector spokeTargetAnglesU = tgt->spokeTargetAnglesU(); + if(verbosityLevel > 0) + std::cout << "Printing information about production target supports:\n"; + + const double targetAngle = tgt->rotHaymanY(); //assume target angle is only in the x-z plane for supports + CLHEP::HepRotation* rodRot = reg.add(CLHEP::HepRotation(CLHEP::HepRotation::IDENTITY)); + rodRot->rotateY(-1.*targetAngle); + + for(int istream = 0; istream < 2; ++istream) { + for(int ispoke = 0; ispoke < nspokesperside; ++ispoke) { + const double wheelAngle = supportWheelRodAngles[ispoke]*CLHEP::degree; + //get angle of the support rod on the wheel and the angle on the target the wire connects to + const double targetWireAngle = (istream == 0) ? spokeTargetAnglesD[ispoke]*CLHEP::degree + : spokeTargetAnglesU[ispoke]*CLHEP::degree; + double rWheel = supportWheelRodRadialOffset[ispoke]; // radius of the wheel to attach to + CLHEP::Hep3Vector rodCenter(rWheel*std::cos(wheelAngle), rWheel*std::sin(wheelAngle), 0.); + const double rodOffset = supportWheelRodOffset[ispoke]; + rodCenter += CLHEP::Hep3Vector(std::sin(targetAngle)*rodOffset, 0., std::cos(targetAngle)*rodOffset); + if(istream == 0) { //only do once + //add the features near the support rods in the bicycle wheel + const double featureAngle = supportWheelFeatureAngles[ispoke]*CLHEP::degree; //angle of feature center + const double featureArc = supportWheelFeatureArcs[ispoke]*CLHEP::degree; //width in angle + const double featureRIn = supportWheelFeatureRIns[ispoke]; //inner radius of feature + const double featureROut = tgt->supportWheelRIn() + smallGap; //ensure they overlap for union + // double featureR = (featureRIn + featureROut)/2.; //radius of feature center + // CLHEP::Hep3Vector featureCenter(featureR*cos(featureAngle), featureR*sin(featureAngle), 0.); + CLHEP::Hep3Vector featureCenter(localWheelCenter); //center is wheel center + double featureParams[] = {featureRIn, featureROut, tgt->supportWheelHL(), featureAngle - featureArc/2. /*phi0*/, featureArc /*dphi*/}; + G4Tubs* featureTubs = new G4Tubs("ProductionTargetSupportFeature_" +std::to_string(ispoke), + featureParams[0], featureParams[1], featureParams[2], featureParams[3], featureParams[4]); + suppWheelInfo.solid = new G4UnionSolid("ProductionTargetSupportWheelFeature_union_"+std::to_string(ispoke), + suppWheelInfo.solid, featureTubs, 0, featureCenter); + //add the support rod to the wheel + G4Tubs* rodTubs = new G4Tubs("ProductionTargetSupportRod_" +std::to_string(ispoke), + 0., supportWheelRodRadius[ispoke], supportWheelRodHL[ispoke], 0., CLHEP::twopi); + suppWheelInfo.solid = new G4UnionSolid("ProductionTargetSupportWheelRod_union_"+std::to_string(ispoke), + suppWheelInfo.solid, rodTubs, rodRot, rodCenter); } - if (verbosityLevel > 0){G4cout << __PRETTY_FUNCTION__ << " ending at z=" << _currentZ << G4endl;} + const int side = (1-2*istream); //+1 or -1 + //info about wire connection + //get end of the rod on this side + CLHEP::Hep3Vector rodAxis(std::sin(targetAngle), 0., std::cos(targetAngle)); + CLHEP::Hep3Vector wheelPos(rodCenter); + wheelPos += side*supportWheelRodHL[ispoke]*rodAxis; + //translate from rod center to edge + CLHEP::Hep3Vector rodCenterToWire(std::cos(wheelAngle)*std::cos(targetAngle), + std::sin(wheelAngle)*std::cos(targetAngle), + -std::cos(wheelAngle)*std::sin(targetAngle)); + wheelPos -= supportWheelRodRadius[ispoke]*rodCenterToWire; + double zWireRodOffset = (istream == 0) ? supportWheelRodWireOffsetD[ispoke] : supportWheelRodWireOffsetU[ispoke]; + wheelPos -= side*zWireRodOffset*rodAxis; + + //get wire position on target + CLHEP::Hep3Vector targetPos(rTarget*cos(targetWireAngle), rTarget*sin(targetWireAngle), side*zTarget); + targetPos = tgt->productionTargetRotation().inverse()*targetPos; //rotate from target frame to mother frame + CLHEP::Hep3Vector spokeAxis((wheelPos-targetPos).unit()); + CLHEP::Hep3Vector targetAxis(0.,0.,side); + targetAxis = tgt->productionTargetRotation().inverse()*targetAxis; + CLHEP::Hep3Vector zax(0.,0.,1.); + if(verbosityLevel > 0) + std::cout << "istream " << istream << " ispoke " << ispoke << std::endl + << "target pos " << targetPos << "\nwheel pos " << wheelPos << std::endl + << "Target axis " << targetAxis << "\nSpoke axis " << spokeAxis << std::endl + << "Rod axis " << rodAxis << "\nRod center to wire axis " << rodCenterToWire << std::endl; + //to remove overlaps where the wire connects, need angle of wire and surface connecting to + //remove overlap at target + double wireTargetAngle = targetAxis.angle(-1.*spokeAxis); + double deltaLength = (std::abs(std::tan(wireTargetAngle)) > 1.e-6) ? std::abs(tgt->spokeRadius()/std::tan(wireTargetAngle)) : 0.; //give up if ~paralle + targetPos += (deltaLength+0.1)*spokeAxis; //subtract off the length + if(verbosityLevel > 0) + std::cout << "wire target angle " << wireTargetAngle << " delta L " << deltaLength + << " target pos " << targetPos <spokeRadius()); + wheelPos -= (deltaLength+1.)*spokeAxis; + if(verbosityLevel > 0) + std::cout << "wire rod angle " << wireRodAngle << " delta L " << deltaLength + << " wheel pos " << wheelPos <spokeRadius(), 0.5*spokeLength); + CLHEP::HepRotation* spokeRot = reg.add(CLHEP::HepRotation(spokeAxis.cross(zax), spokeAxis.angle(zax))); + std::stringstream spokeName; + spokeName << "ProductionTargetSpokeWire_" ; + if(istream == 0) + spokeName << "Downstream_"; + else + spokeName << "Upstream_"; + spokeName << ispoke; + + VolumeInfo spokeInfo = nestTubs( spokeName.str(), + spokeParams, + spokeMaterial, + spokeRot, + spokeCenter, + prodTargetMotherInfo, + 0, + G4Colour::Gray(), + "PS" + ); + + } //end spokes loop + } //end stream loop + finishNesting(suppWheelInfo, + suppWheelMaterial, + 0, + localWheelCenter, + prodTargetMotherInfo.logical, + 0, + G4Colour::Gray(), + "PS" + ); + + } //end adding support structures + } //end constructHaymanTarget + + void constructStickmanTarget(VolumeInfo const & parent, SimpleConfig const & _config) { + Mu2eG4Helper & _helper = *(art::ServiceHandle()); + AntiLeakRegistry & reg = _helper.antiLeakRegistry(); + + int verbosityLevel = _config.getInt("PSStickman.verbosityLevel"); + verbosityLevel >0 && + cout << __func__ << " verbosityLevel on Stickman1.0 : " << verbosityLevel << endl; + G4ThreeVector _hallOriginInMu2e = parent.centerInMu2e(); + // Create variable to avoid multiple look-up + + G4GeometryOptions* geomOptions = art::ServiceHandle()->geomOptions(); + + bool prodTargetVisible = geomOptions->isVisible( "ProductionTarget" ); + bool prodTargetSolid = geomOptions->isSolid ( "ProductionTarget" ); + bool forceAuxEdgeVisible = geomOptions->forceAuxEdgeVisible( "ProductionTarget" ); + bool placePV = geomOptions->placePV( "ProductionTarget" ); + bool doSurfaceCheck = geomOptions->doSurfaceCheck( "ProductionTarget" ); + + // begin all names with ProductionTarget so when we build sensitive detectors we can make them all + // sensitive at once with LVname.find("ProductionTarget") !=std::string::npos + // + // Build the production target. + + GeomHandle tgt; + // allows a vector of plate materials so each plate can be different. Initialize to have same number of + // entries as the number of plates. Fin and core materials are the same for each plate. + vector prodTargetPlateMaterials(tgt->numberOfPlates(), nullptr); + for (int i = 0; i < tgt->numberOfPlates(); ++i) { + prodTargetPlateMaterials.at(i) = findMaterialOrThrow(tgt->plateMaterial(i)); + } + G4Material* prodTargetRodMaterial = findMaterialOrThrow(tgt->rodMaterial()); + G4Material* prodTargetSpacerMaterial = findMaterialOrThrow(tgt->spacerMaterial()); + G4Material* prodTargetSupportRingMaterial = findMaterialOrThrow(tgt->stickmanSupportRingMaterial()); + + TubsParams prodTargetMotherParams( 0. + ,tgt->productionTargetMotherOuterRadius() + ,tgt->productionTargetMotherHalfLength()); + + G4ThreeVector _localCenter(0.0,0.0,0.0); + G4RotationMatrix* targetRotation = reg.add(G4RotationMatrix(tgt->productionTargetRotation().inverse())); + if (verbosityLevel > 2){G4cout << __PRETTY_FUNCTION__ << "target rotation = " << *targetRotation << G4endl;} + VolumeInfo prodTargetMotherInfo = nestTubs( "ProductionTargetMother", + prodTargetMotherParams, + parent.logical->GetMaterial(), + 0, + tgt->stickmanProdTargetPosition() - parent.centerInMu2e(), + parent, + 0, + G4Colour::Blue(), + "PS" + ); + + if (verbosityLevel > 0){ + G4cout << "target position and hall origin = " + << tgt->stickmanProdTargetPosition() << "\n" + << _hallOriginInMu2e << " " << parent.centerInMu2e() << G4endl; + } + if (verbosityLevel > 2){ + G4cout << __PRETTY_FUNCTION__ << "created prodTargetMotherInfo " + << tgt->productionTargetMotherOuterRadius() + << " " <productionTargetMotherHalfLength() << G4endl; + } + + // tiny offset to avoid precision issues + constexpr double stickmanMagicOffset = 0.0001; + + // store the fin rotations for later use + // since used for G4UnionSolid construction, a passive rotation is needed, hence the negative sign on the angle. + // These are the rotations to go from the plate frame to the fin frame + std::vector stickmanFinRotations; + for (int ithFin = 0; ithFin < tgt->nStickmanFins(); ++ithFin) { + CLHEP::HepRotation* finRotation = reg.add(CLHEP::HepRotation::IDENTITY); + finRotation->rotateZ(-tgt->plateFinAngle(ithFin)); + stickmanFinRotations.emplace_back(finRotation); + } + + // cache the plate solids + std::map stickmanPlateSolidCache; + // lambda to get the plate solid for a given plate. The solid is cached based on the plate parameters that affect the shape. + // The cache key encodes these parameters; if a solid is not found in the cache, it is created and added to the cache before + //being returned. + // Here chose to construct plates as a single union solid with different materials for each plate, + // rather than separate fins, cores, and lugs, because the fins is much wider than the Hayman design. + // If separated, missing material between fins and cores would be large. + // Note that the geometry center of a G4UnionSolid is the same as the first solid in the union, + // and I start with the plate core and then union on the fins and lugs. This means that when placing the + // plates the position will be based on the center of the core. + auto getStickmanPlateSolid = [&](int ithPlate) -> G4VSolid* { + const double plateFilletRadius = + (tgt->addFilletToPlateCore() || tgt->addFilletToPlateLug()) ? tgt->plateFilletRadius() : 0.; + const std::string cacheKey = std::string("StickmanPlate_") + + std::to_string(tgt->plateROut(ithPlate)) + "_" // Core radius + + std::to_string(tgt->plateThickness(ithPlate)) + "_" // plate thickness + + std::to_string(tgt->plateLugThickness(ithPlate)) + "_" // lug thickness + + std::to_string(static_cast(tgt->addFilletToPlateCore())) + "_" // whether to add fillet to plate core + + std::to_string(static_cast(tgt->addFilletToPlateLug())) + "_" // whether to add fillet to plate lug + + std::to_string(plateFilletRadius); // fillet radius (if applicable) + + // check if the solid is already in the cache + auto cached = stickmanPlateSolidCache.find(cacheKey); + if (cached != stickmanPlateSolidCache.end()) { + if (verbosityLevel > 3) { + G4cout << __PRETTY_FUNCTION__ << " reusing cached plate solid for plate " + << ithPlate << G4endl; } - ++coreCopyNumber; - // - // if this is the last section, add on one more beginning block. decided not to extend vectors and have zero segments, just ugly - if (ithSection == (numberOfSections -1) ){ - startName = "ProductionTargetStartingCoreSection_" + std::to_string(ithSection+1+1); - TubsParams startingSegmentParamsEnd(0.,targetRadius,tgt->startingSectionThickness().at(ithSection)/2.); + return cached->second; + } - // - //set currentZ to be in the center of the starting section - _currentZ += tgt->startingSectionThickness().at(ithSection)/2.; - if (verbosityLevel > 0){G4cout << __PRETTY_FUNCTION__ << " ithSection+1 , startingSegmentCenter at " << ithSection+1 << " " <<_currentZ << G4endl;} - _startingSegmentCenter.setZ(_currentZ); - CLHEP::Hep3Vector currentStartingSegmentCenter = tgt->productionTargetRotation().inverse()*_startingSegmentCenter; - double currentHalfStartingSegment = tgt->startingSectionThickness().at(ithSection)/2.; - startingSegment = nestTubs(startName, - startingSegmentParamsEnd, - prodTargetCoreMaterial, - &tgt->productionTargetRotation(), - currentStartingSegmentCenter, - prodTargetMotherInfo, - ithSection+1, - prodTargetVisible, - G4Colour::Yellow(), - prodTargetSolid, - forceAuxEdgeVisible, - placePV, - doSurfaceCheck - ); - - //build fins around starting segment - for (int ithFin = 0; ithFin < tgt->nHaymanFins(); ++ithFin){ - - double currentFinAngle = tgt->finAngles().at(ithFin); - const std::string name = "ProductionTargetFinStartingSection_" + std::to_string(ithSection+1+1) + "_Fin_" + std::to_string(ithFin); + // if not in cache, create the solid and add to cache + const double plateHalfThickness = tgt->plateThickness(ithPlate)/2.; + const double lugHalfThickness = tgt->plateLugThickness(ithPlate)/2.; + const double finHalfWidth = tgt->plateFinWidth()/2.; + const double finHalfLength = tgt->plateFinOuterRadius()/2.; + const double lugZShift = lugHalfThickness - plateHalfThickness; + if (verbosityLevel > 3) { + G4cout << __PRETTY_FUNCTION__ << " creating new plate solid for plate " + << ithPlate << " with cache key " << cacheKey << G4endl; + } - if (verbosityLevel > 1) - {G4cout << "Fin Section and Angle = " << name << " " << currentFinAngle << " fin copy number = " << finCopyNumber << G4endl;} + // solid for the plate core + G4VSolid* plateSolid = reg.add(new G4Tubs(cacheKey + "_Core", // + 0., // inner radius + tgt->plateROut(ithPlate), // outer radius + plateHalfThickness, // z half length + 0., // starting angle + CLHEP::twopi)); // segment angle is full circle + if (verbosityLevel > 3) { + G4cout << __PRETTY_FUNCTION__ << " core solid parameters (rOut, halfLength) = (" + << tgt->plateROut(ithPlate) << ", " << plateHalfThickness << ")" << G4endl; + } - double finHalfHeightAboveTarget = (tgt->finOuterRadius() - tgt->rOut())/2.; - G4Box* finBox = new G4Box("finBox",tgt->haymanFinThickness()/2.,finHalfHeightAboveTarget,currentHalfStartingSegment); - // - // this tubs is in the xy plane with an extent along z. So phi goes in the xy plane. - // need extra length for visualization tool and overlap avoidance - G4Tubs* tubCutout = new G4Tubs("finCutout",0.,tgt->rOut(),currentHalfStartingSegment + magicOffset,dSphi,dPphi); + // solid for the fin + G4Box* finBox = reg.add(new G4Box(cacheKey + "_FinBox", // + finHalfLength, // x half length + finHalfWidth, // y half length + plateHalfThickness)); // z half length + if (verbosityLevel > 3) { + G4cout << __PRETTY_FUNCTION__ << " fin solid parameters (halfLength, halfWidth, halfThickness) = (" + << finHalfLength << ", " << finHalfWidth << ", " << plateHalfThickness << ")" << G4endl; + } - CLHEP::HepRotation* rotFin = reg.add(CLHEP::HepRotation(tgt->productionTargetRotation().inverse())); - CLHEP::HepRotation* rotFinG4 = reg.add(CLHEP::HepRotation(tgt->productionTargetRotation())); - // - // ok here I have a little confusing fix. The fin is built along the y-axis. But the rotation is given wrt the x axis. Hence I need to - // subtract off that 90^o - double rotAdj = -M_PI/2 + currentFinAngle; - rotFin->rotateZ(rotAdj); - // - // g4 version. - rotFinG4->rotateZ(-rotAdj); - ++finCopyNumber; + // solid for the lug + G4Tubs* lugSolid = reg.add(new G4Tubs(cacheKey + "_Lug", // + tgt->plateLugInnerRadius(), // inner radius + tgt->plateLugOuterRadius(), // outer radius + lugHalfThickness, // z half length + 0., // starting angle + CLHEP::twopi)); // segment angle is full circle + if (verbosityLevel > 3) { + G4cout << __PRETTY_FUNCTION__ << " lug solid parameters (rIn, rOut, halfThickness) = (" + << tgt->plateLugInnerRadius() << ", " << tgt->plateLugOuterRadius() << ", " << lugHalfThickness << ")" << G4endl; + } - // - // for debugging - CLHEP::Hep3Vector offset(0.,0.,0.); - // - // rotate the offset vector - G4ThreeVector finTranslation = currentStartingSegmentCenter; - // - // G4 translates and then rotates. I want to offset the fin to its final location relative to the core and then let G4 perform rotations, - // both of the fin about the z axis and then the entire production target rotation angle - double distanceFromCenter = finHalfHeightAboveTarget + tgt->rOut(); - CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*cos(currentFinAngle),distanceFromCenter*sin(currentFinAngle),0.); - // - // finTranslation was already put in the target rotated frame. I have to rotate the offset vector as well - offsetRelativeToCore *= tgt->productionTargetRotation().inverse(); - finTranslation = finTranslation + offsetRelativeToCore; - if (verbosityLevel > 2){G4cout << "finTranslation = " << finTranslation << G4endl;} - CLHEP::Hep3Vector downshift = CLHEP::Hep3Vector(0.,-finHalfHeightAboveTarget - targetRadius*cos(angularSize), 0.); - - VolumeInfo finWithCutoutInfo(name,finTranslation,prodTargetMotherInfo.centerInWorld); - finWithCutoutInfo.solid = new G4SubtractionSolid(name,finBox,tubCutout,nullptr,downshift); - finishNesting(finWithCutoutInfo - ,prodTargetFinMaterial - ,rotFinsG4.at(ithFin) - ,finTranslation - ,prodTargetMotherInfo.logical - ,finCopyNumber - ,G4Colour::White() - ,"ProductionTarget" - ,verbosityLevel>1); + // add fillets to the plate core + // for a fin at 0 degree start with unioning a G4ExtrudedSolid with a triangular cross section that has vertices in + // the x-y plane at (0,0), + // (sqrt((rOut+FilletRadius)^2-(plateFinWidth/2+FilletRadius)^2), plateFinWidth/2+FilletRadius), + // (sqrt((rOut+FilletRadius)^2-(plateFinWidth/2+FilletRadius)^2), -plateFinWidth/2-FilletRadius), + // and then extrude along z by plateThickness(ithPlate) starting at z = _currentZ. The extrusion has no twist and + // a scale of 1.0 at both ends. + // The fillet is then completed by subtracting a cylinder with radius of the fillet radius (plus a small tolerance) + // and height plateThickness(ithPlate) that is centered at + // (sqrt((rOut+FilletRadius)^2-(plateFinWidth/2+FilletRadius)^2), plateFinWidth/2+FilletRadius), + // and the other fillet is completed by subtracting a similar cylinder at + // (sqrt((rOut+FilletRadius)^2-(plateFinWidth/2+FilletRadius)^2), -plateFinWidth/2-FilletRadius), + // For fins at other angles, the same is done but the vertices are rotated by the fin angle and the centers of the + // cylinders for subtraction are also rotated by the fin angle. + G4VSolid* plateCoreFillet = nullptr; + if (tgt->addFilletToPlateCore()) { + // get the coordinates of the fillet union vertices for the fin at 0 degree + const double coreFilletX = std::sqrt(std::pow(tgt->plateROut(ithPlate) + plateFilletRadius, 2) + - std::pow(finHalfWidth + plateFilletRadius, 2)); + std::vector coreFilletVertices = { + G4TwoVector(0., 0.), + G4TwoVector(coreFilletX, finHalfWidth + plateFilletRadius), + G4TwoVector(coreFilletX, -finHalfWidth - plateFilletRadius) + }; + // create the extrusion solid for the fillet + G4ExtrudedSolid* coreFilletExtrusion = reg.add(new G4ExtrudedSolid( + cacheKey + "_CoreFilletExtrusion", // + coreFilletVertices, // vertices + plateHalfThickness, // half length in z + G4TwoVector(0.,0.), // shift of polygon center at -half length + 1.0, // scale of polygon at -half length + G4TwoVector(0.,0.), // shift of polygon center at +half length + 1.0)); // scale of polygon at +half length + if (verbosityLevel > 3) { + G4cout << __PRETTY_FUNCTION__ << " core fillet extrusion vertices = (" + << coreFilletVertices.at(0) << ", " + << coreFilletVertices.at(1) << ", " + << coreFilletVertices.at(2) << ")" << G4endl; + G4cout << __PRETTY_FUNCTION__ << " core fillet extrusion half length = " << plateHalfThickness << G4endl; + } - } - // - // set current z to be at the end of this starting segment as a final check - _currentZ += tgt->startingSectionThickness().at(ithSection)/2.; - if (verbosityLevel > 0){ - G4cout << __PRETTY_FUNCTION__ << " ithSection+1 , startingSegment Z ends at " << ithSection+1 << " " <<_currentZ << G4endl; - G4cout << " with copy number for last segment = " << ithSection+1 << G4endl; - } - /***********and now the final end ring ***********/ - std::string name = "ProductionTargetPositiveEndRing"; - G4Tubs* supportRing = new G4Tubs(name, - tgt->supportRingInnerRadius(), - tgt->supportRingOuterRadius(), - tgt->supportRingLength()/2., - 0., - 2.*M_PI); - if (verbosityLevel > 0){ - G4cout << __PRETTY_FUNCTION__ << ": \n " << name.c_str() << ":\n (rin, rout, half length) = (" - << tgt->supportRingInnerRadius() << ", " << tgt->supportRingOuterRadius() - << ", " << tgt->supportRingLength()/2. << ")" << G4endl; - } - // move ring to end of target and then back by cutout size (signs for positive side) - G4ThreeVector ringTranslation = currentStartingSegmentCenter - + tgt->productionTargetRotation().inverse()*CLHEP::Hep3Vector(0.,0.,tgt->supportRingLength()/2.); + // create the cutout cylinder for the fillet + G4Tubs* coreFilletCutout = reg.add(new G4Tubs( + cacheKey + "_CoreFilletCutout", // + 0., // inner radius + plateFilletRadius + stickmanMagicOffset, // outer radius (add small offset to avoid precision issues) + plateHalfThickness + stickmanMagicOffset, // half length in z + 0., // starting angle + CLHEP::twopi)); // segment angle is full circle + // make 2 cutouts for the fillet, one at the upper edge of the fin and one at the lower edge of the fin + G4VSolid* coreFilletUpper = reg.add(new G4SubtractionSolid( + cacheKey + "_CoreFilletUpper", // + coreFilletExtrusion, // base solid + coreFilletCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(coreFilletX, + finHalfWidth + plateFilletRadius, + 0.))); // translation of cutout solid + plateCoreFillet = reg.add(new G4SubtractionSolid( + cacheKey + "_CoreFilletLower", // + coreFilletUpper, // base solid + coreFilletCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(coreFilletX, + -finHalfWidth - plateFilletRadius, + 0.))); // translation of cutout solid + if (verbosityLevel > 3) { + G4cout << __PRETTY_FUNCTION__ << " core fillet cutout radius = " << plateFilletRadius + stickmanMagicOffset << G4endl; + G4cout << __PRETTY_FUNCTION__ << " core fillet upper cutout center = (" << coreFilletX << ", " << finHalfWidth + plateFilletRadius << ", 0)" << G4endl; + G4cout << __PRETTY_FUNCTION__ << " core fillet lower cutout center = (" << coreFilletX << ", " << -finHalfWidth - plateFilletRadius << ", 0)" << G4endl; + } + } // end of plate core fillet construction + + // add fillets to the plate lugs + // for a fin at 0 degree, start with unioning a G4ExtrudedSolid with a triangular cross section that has vertices in + // the x-y plane at (plateCenterToLugCenter, 0), + // (plateCenterToLugCenter-sqrt((plateLugOuterRadius+FilletRadius)^2-(plateFinWidth/2+FilletRadius)^2), plateFinWidth/2+FilletRadius), + // (plateCenterToLugCenter-sqrt((plateLugOuterRadius+FilletRadius)^2-(plateFinWidth/2+FilletRadius)^2), -plateFinWidth/2-FilletRadius), + // and then extrude along z by plateThickness(ithPlate) starting at z = _currentZ. The extrusion has no twist and + // a scale of 1.0 at both ends. + + // The fillet is then completed by subtracting a cylinder with radius of the fillet radius (plus a small tolerance) + // and height plateThickness(ithPlate) that is centered at + // (plateCenterToLugCenter-sqrt((plateLugOuterRadius+FilletRadius)^2-(plateFinWidth/2+FilletRadius)^2), plateFinWidth/2+FilletRadius), + // and the other fillet is completed by subtracting a similar cylinder at + // (plateCenterToLugCenter-sqrt((plateLugOuterRadius+FilletRadius)^2-(plateFinWidth/2+FilletRadius)^2), -plateFinWidth/2-FilletRadius), + // For fins at other angles, the same is done but the vertices are rotated by the fin angle and the centers of the + // cylinders for subtraction are also rotated by the fin angle. + G4VSolid* plateLugFillet = nullptr; + if (tgt->addFilletToPlateLug()) { + // get the coordinates of the fillet union vertices for the fin at 0 degree + const double lugFilletX = tgt->plateCenterToLugCenter() + - std::sqrt(std::pow(tgt->plateLugOuterRadius() + plateFilletRadius, 2) + - std::pow(finHalfWidth + plateFilletRadius, 2)); + std::vector lugFilletVertices = { + G4TwoVector(tgt->plateCenterToLugCenter(), 0.), + G4TwoVector(lugFilletX, finHalfWidth + plateFilletRadius), + G4TwoVector(lugFilletX, -finHalfWidth - plateFilletRadius) + }; + // create the extrusion solid for the fillet + G4ExtrudedSolid* lugFilletExtrusion = reg.add(new G4ExtrudedSolid( + cacheKey + "_LugFilletExtrusion", // + lugFilletVertices, // vertices + plateHalfThickness, // half length in z + G4TwoVector(0.,0.), // shift of polygon center at -half length + 1.0, // scale of polygon at -half length + G4TwoVector(0.,0.), // shift of polygon center at +half length + 1.0)); // scale of polygon at +half length + if (verbosityLevel > 3) { + G4cout << __PRETTY_FUNCTION__ << " lug fillet extrusion vertices = (" + << lugFilletVertices.at(0) << ", " + << lugFilletVertices.at(1) << ", " + << lugFilletVertices.at(2) << ")" << G4endl; + G4cout << __PRETTY_FUNCTION__ << " lug fillet extrusion half length = " << plateHalfThickness << G4endl; + } - if (verbosityLevel > 0){ - G4cout << " center = (0., 0., " << (currentStartingSegmentCenter.mag() + tgt->supportRingLength()/2.) << ")" << G4endl; - G4cout << " center (wrt target mother) = " << ringTranslation << G4endl; - } - std::vector ringWithCutoutSolid; - std::string nameRing = "ProductionTargetPositiveRingCutout"; + // create the cutout cylinder for the fillet + G4Tubs* lugFilletCutout = reg.add(new G4Tubs( + cacheKey + "_LugFilletCutout", // + 0., // inner radius + plateFilletRadius + stickmanMagicOffset, // outer radius (add small offset to avoid precision issues) + plateHalfThickness + stickmanMagicOffset, // half length in z + 0., // starting angle + CLHEP::twopi)); // segment angle is full circle + // lug center needs to be cut out too. Create a second cutout cylinder for the lug center with radius of + // plateLugInnerRadius and height of plateThickness(ithPlate) centered at (plateCenterToLugCenter, 0) + G4Tubs* lugCenterCutout = reg.add(new G4Tubs( + cacheKey + "_LugCenterCutout", // + 0., // inner radius + tgt->plateLugInnerRadius() + stickmanMagicOffset, // outer radius (add small offset to avoid precision issues) + plateHalfThickness + stickmanMagicOffset, // half length in z + 0., // starting angle + CLHEP::twopi)); // segment angle is full circle + // cut out the lug center from the extrusion first + G4VSolid* lugFilletWithCenterCutout = reg.add(new G4SubtractionSolid( + cacheKey + "_LugFilletWithCenterCutout", // + lugFilletExtrusion, // base solid + lugCenterCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(tgt->plateCenterToLugCenter(), 0., 0.))); // translation of cutout solid + // make 2 cutouts for the fillet, one at the upper edge of the fin and one at the lower edge of the fin + G4VSolid* lugFilletUpper = reg.add(new G4SubtractionSolid( + cacheKey + "_LugFilletUpper", // + lugFilletWithCenterCutout, // base solid + lugFilletCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(lugFilletX, + finHalfWidth + plateFilletRadius, + 0.))); // translation of cutout solid + plateLugFillet = reg.add(new G4SubtractionSolid( + cacheKey + "_LugFilletLower", // + lugFilletUpper, // base solid + lugFilletCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(lugFilletX, + -finHalfWidth - plateFilletRadius, + 0.))); // translation of cutout solid + if (verbosityLevel > 3) { + G4cout << __PRETTY_FUNCTION__ << " lug fillet cutout radius = " << plateFilletRadius + stickmanMagicOffset << G4endl; + G4cout << __PRETTY_FUNCTION__ << " lug center cutout radius = " << tgt->plateLugInnerRadius() + stickmanMagicOffset << G4endl; + G4cout << __PRETTY_FUNCTION__ << " lug fillet upper cutout center = (" << lugFilletX << ", " << finHalfWidth + plateFilletRadius << ", 0)" << G4endl; + G4cout << __PRETTY_FUNCTION__ << " lug fillet lower cutout center = (" << lugFilletX << ", " << -finHalfWidth - plateFilletRadius << ", 0)" << G4endl; + } + } // end of plate lug fillet construction - for (int ithFin = 0; ithFin < tgt->nHaymanFins(); ++ithFin){ - double currentFinAngle = tgt->finAngles().at(ithFin); - const std::string name = "ProductionTargetPositiveRingCutout_" + std::to_string(ithSection) - + "_Segment_" +"_Fin_" + std::to_string(ithFin); - ++boxCopyNumber; + // + // assemble the plate by unioning the core, fins, and lugs (with fillets if applicable) + for (int ithFin = 0; ithFin < tgt->nStickmanFins(); ++ithFin) { + const double currentFinAngle = tgt->plateFinAngle(ithFin); + // fin and lug shifts of the current fin + const CLHEP::Hep3Vector finShift(finHalfLength * std::cos(currentFinAngle), + finHalfLength * std::sin(currentFinAngle), + 0.); + const CLHEP::Hep3Vector lugShift(tgt->plateCenterToLugCenter() * std::cos(currentFinAngle), + tgt->plateCenterToLugCenter() * std::sin(currentFinAngle), + lugZShift); + // note the lug is shifted in z to be flush with the edge of the plate core + + // add fin + plateSolid = reg.add(new G4UnionSolid(cacheKey + "_FinUnion_" + std::to_string(ithFin), + plateSolid, // base solid + finBox, // union solid + stickmanFinRotations.at(ithFin), // rotation of union solid + finShift)); // translation of union solid + if (verbosityLevel > 3) { + G4cout << __PRETTY_FUNCTION__ << " fin " << ithFin << " added with angle of " << tgt->plateFinAngle(ithFin)/CLHEP::degree << " deg" << G4endl; + } + // add core fillets if applicable + if (plateCoreFillet != nullptr) { + plateSolid = reg.add(new G4UnionSolid(cacheKey + "_CoreFilletUnion_" + std::to_string(ithFin), + plateSolid, // base solid + plateCoreFillet, // union solid + stickmanFinRotations.at(ithFin), // rotation of union solid + CLHEP::Hep3Vector(0.,0.,0.))); // translation of union solid + } - // - // G4 translates and then rotates. I want to offset the fin to its final location relative to the core and then let G4 perform rotations, - // both of the fin about the z axis and then the entire production target rotation angle - double distanceFromCenter = cutoutHeightAboveTarget; - - CLHEP::Hep3Vector boxShift(distanceFromCenter*cos(currentFinAngle),distanceFromCenter*sin(currentFinAngle),0.); - CLHEP::Hep3Vector offsetRelativeToCore(distanceFromCenter*cos(currentFinAngle),distanceFromCenter*sin(currentFinAngle), - -tgt->supportRingCutoutLength()/2. - magicOffset);// note flipped sign from other end of target - - CLHEP::Hep3Vector downshift = CLHEP::Hep3Vector(0.,0.,0.);//for debugging - offsetRelativeToCore = offsetRelativeToCore + downshift; - if (ithFin == 0){ - ringWithCutoutSolid.push_back(new G4SubtractionSolid(name,supportRing,cutoutBox,rotBoxesG4.at(ithFin),offsetRelativeToCore)); - }else { - ringWithCutoutSolid.push_back(new G4SubtractionSolid(name,ringWithCutoutSolid.at(ithFin-1),cutoutBox,rotBoxesG4.at(ithFin),offsetRelativeToCore)); - } - } + // add lug + plateSolid = reg.add(new G4UnionSolid(cacheKey + "_LugUnion_" + std::to_string(ithFin), + plateSolid, // base solid + lugSolid, // union solid + nullptr, // rotation of union solid + lugShift)); // translation of union solid + if (verbosityLevel > 3) { + G4cout << __PRETTY_FUNCTION__ << " lug " << ithFin << " added with shift (" << lugShift.x() << ", " << lugShift.y() << ", " << lugShift.z() << ")" << G4endl; + } + // add lug fillets if applicable + if (plateLugFillet != nullptr) { + plateSolid = reg.add(new G4UnionSolid(cacheKey + "_LugFilletUnion_" + std::to_string(ithFin), + plateSolid, // base solid + plateLugFillet, // union solid + stickmanFinRotations.at(ithFin), // rotation of union solid + CLHEP::Hep3Vector(0.,0.,0.))); // translation of union solid + } + } // end of loop over fins to assemble the plate solid - if (verbosityLevel > 0){ - G4cout << " center (at end) = (0., 0., " << ringTranslation.mag() << ")" <plateThickness(ithPlate)/2.; + const double plateCenterZ = _currentZ + plateHalfThickness; // this is the z position of the core + const CLHEP::Hep3Vector localPlateCenter(0.,0.,plateCenterZ); // this is in the target frame + const CLHEP::Hep3Vector currentPlateCenter = tgt->productionTargetRotation().inverse()*localPlateCenter; + // productionTargetRotation() rotates to the target frame, so to get the position in the world frame we + // apply the inverse rotation to the local plate center. + + if (verbosityLevel > 0) { + G4cout << __PRETTY_FUNCTION__ + << " plate " << ithPlate + << " startZ=" << _currentZ + << " centerZ=" << plateCenterZ + << " endZ=" << (_currentZ + tgt->plateLugThickness(ithPlate)) + << " material=" << tgt->plateMaterial(ithPlate) + << G4endl; + } + if (verbosityLevel > 2) { + G4cout << __PRETTY_FUNCTION__ + << " rOut=" << tgt->plateROut(ithPlate) + << " thickness=" << tgt->plateThickness(ithPlate) + << " finWidth=" << tgt->plateFinWidth() + << " finReach=" << tgt->plateFinOuterRadius() + << " lugThickness=" << tgt->plateLugThickness(ithPlate) + << " lug(inner,outer)=(" + << tgt->plateLugInnerRadius() << "," + << tgt->plateLugOuterRadius() << ")" + << " addFilletToPlateCore=" << tgt->addFilletToPlateCore() + << " addFilletToPlateLug=" << tgt->addFilletToPlateLug() + << G4endl; + } + + // get the solid for the plate, using the lambda defined above + G4VSolid* plateSolid = getStickmanPlateSolid(ithPlate); + + // add the plate to the target mother + VolumeInfo plateInfo(plateName, currentPlateCenter, prodTargetMotherInfo.centerInWorld); + plateInfo.solid = plateSolid; + finishNesting(plateInfo // volume info as defined above + ,prodTargetPlateMaterials.at(ithPlate) // material for the plate, can be different for each plate + ,&tgt->productionTargetRotation() // PASSIVE rotation + ,currentPlateCenter // translation to place the plate at the correct z position within the target mother + ,prodTargetMotherInfo.logical // parent volume logical to place the plate in + ,ithPlate // copy number + ,G4Colour::Yellow() // color for visualization + ,"ProductionTarget" // lookup token + ,verbosityLevel>1); // pass the verbosity level for printing info in finishNesting + + // update _currentZ to the END of the next plate, which is the start of the current plate + plateLugThickness(ithPlate) + _currentZ += tgt->plateLugThickness(ithPlate); + } // end of loop over plates + + // add the rods and spacers + G4VSolid* spacerSolid = reg.add(new G4Tubs( + "ProductionTargetSpacerSolid", // + tgt->spacerInnerRadius(), // inner radius + tgt->spacerOuterRadius(), // outer radius + tgt->spacerHalfLength(), // half length + 0., // start angle + CLHEP::twopi)); // sweep angle + const double rodHalfLength = tgt->rodHalfLength() - stickmanMagicOffset; // make it slightly shorter to avoid precesion issues + const double spacerCenterZOffset = tgt->rodHalfLength() - tgt->spacerHalfLength(); + + for (int ithRod = 0; ithRod < tgt->nStickmanFins(); ++ithRod) { + const double rodAngle = tgt->plateFinAngle(ithRod); + const CLHEP::Hep3Vector rodCenterInTargetFrame( + tgt->plateCenterToLugCenter() * std::cos(rodAngle), + tgt->plateCenterToLugCenter() * std::sin(rodAngle), + 0.); + const CLHEP::Hep3Vector rodCenter = tgt->productionTargetRotation().inverse() * rodCenterInTargetFrame; + + TubsParams rodParams(0., // rIn + tgt->rodRadius(), // rOut + rodHalfLength); // half length + VolumeInfo rodInfo = nestTubs("ProductionTargetRod_" + std::to_string(ithRod), + rodParams, + prodTargetRodMaterial, + &tgt->productionTargetRotation(), // passive rotation + rodCenter, // translation + prodTargetMotherInfo, // parent volume info for the target mother + ithRod, // copy number + prodTargetVisible, // visibility + G4Colour::Green(), // color for visualization + prodTargetSolid, + forceAuxEdgeVisible, + placePV, + doSurfaceCheck + ); + if (verbosityLevel > 0) { + G4cout << __PRETTY_FUNCTION__ << " rod " << ithRod + << " center=" << rodCenter + << " material=" << tgt->rodMaterial() + << G4endl; + } + + // place the spacers one spacerHalfLength inward from each rod end + const CLHEP::Hep3Vector spacerCenterNegZInTargetFrame( + rodCenterInTargetFrame.x(), + rodCenterInTargetFrame.y(), + -spacerCenterZOffset); + const CLHEP::Hep3Vector spacerCenterPosZInTargetFrame( + rodCenterInTargetFrame.x(), + rodCenterInTargetFrame.y(), + spacerCenterZOffset); + // transform the spacer centers to the parent frame + const CLHEP::Hep3Vector spacerShiftNegZ = + tgt->productionTargetRotation().inverse() * spacerCenterNegZInTargetFrame; + const CLHEP::Hep3Vector spacerShiftPosZ = + tgt->productionTargetRotation().inverse() * spacerCenterPosZInTargetFrame; + if (verbosityLevel > 2) { + G4cout << __PRETTY_FUNCTION__ << " spacer " << ithRod + << " target-frame centers neg/pos=" + << spacerCenterNegZInTargetFrame << " / " + << spacerCenterPosZInTargetFrame + << " mother-frame centers neg/pos=" + << spacerShiftNegZ << " / " + << spacerShiftPosZ + << G4endl; + } + + // add the spacers + std::string spacerNameNegZ = "ProductionTargetSpacerNegZ_" + std::to_string(ithRod); + VolumeInfo spacerInfoNegZ(spacerNameNegZ, spacerShiftNegZ, prodTargetMotherInfo.centerInWorld); + spacerInfoNegZ.solid = spacerSolid; + finishNesting(spacerInfoNegZ, // volume info as defined above + prodTargetSpacerMaterial, // material for the spacer + &tgt->productionTargetRotation(), // rotation + spacerShiftNegZ, // translation + prodTargetMotherInfo.logical, // parent logical volume + 2*ithRod, // copy number + G4Colour::Cyan(), // color for visualization + "ProductionTarget", // lookup token + verbosityLevel>1); // pass the verbosity level for printing info in finishNesting + std::string spacerNamePosZ = "ProductionTargetSpacerPosZ_" + std::to_string(ithRod); + VolumeInfo spacerInfoPosZ(spacerNamePosZ, spacerShiftPosZ, prodTargetMotherInfo.centerInWorld); + spacerInfoPosZ.solid = spacerSolid; + finishNesting(spacerInfoPosZ, // volume info as defined above + prodTargetSpacerMaterial, // material for the spacer + &tgt->productionTargetRotation(), // rotation + spacerShiftPosZ, // translation + prodTargetMotherInfo.logical, // parent logical volume + 2*ithRod + 1, // copy number + G4Colour::Cyan(), // color for visualization + "ProductionTarget", // lookup token + verbosityLevel>1); // pass the verbosity level for printing info in finishNesting + if (verbosityLevel > 0) { + G4cout << __PRETTY_FUNCTION__ << " spacer for rod " << ithRod + << " negZ center=" << spacerShiftNegZ + << " posZ center=" << spacerShiftPosZ + << " material=" << tgt->spacerMaterial() + << G4endl; + } + } // end of loop over rods and spacers + + // now add the two end rings + + // then two cutouts are made by subtracting cylinders with radius of the fillet radius (plus a small tolerance) and height of stickmanSupportRingLength (plus a small tolerance) + // centered at (sqrt((supportRingOuterRadius+filletRadius)^2-(supportRingLugOuterRadius+filletRadius)^2), supportRingLugOuterRadius+filletRadius), and + // (sqrt((supportRingOuterRadius+filletRadius)^2-(supportRingLugOuterRadius+filletRadius)^2), -supportRingLugOuterRadius-filletRadius) + // an additional cutout is made for the inner part of the ring by subtracting a tube with inner radius of 0, outer radius of supportRingInnerRadius, and height of stickmanSupportRingLength (plus a small tolerance) centered at (0,0). + // if addCutoutToSupportRing is true, then additional circular cutouts (through holes) are made. + // The number of cutouts is given by nSupportRingCutouts at angles defined by supportRingCutoutAngles, and the cutouts are made by subtracting cylinders with radius of + // supportRingCutoutInnerRadius. The centers of the cutouts are at supportRingCutoutOffset from the interface between the spacers and the end ring, and at the angles defined + // by supportRingCutoutAngles. The angle is between the through hole and the local radius, the through hole points outward and towards the middle. + + // lambda to build the end support ring solid with lugs and fillets if applicable + // Note the end plate has a flat face on the inside that is flush with the end of the rods / spacers, details like chamfers at + // the rod sockets and those at the cutouts (through holes) are ignored for simplicity. + auto buildStickmanSupportRingSolid = [&](const std::string& tag, bool positiveEnd) -> G4VSolid* { + const double ringHalfLength = tgt->stickmanSupportRingLength()/2.; + // dimensions for the lug stem + const double lugStemHalfLength = 0.5*(tgt->plateCenterToLugCenter() - tgt->stickmanSupportRingInnerRadius()); + const double lugStemCenterRadius = tgt->stickmanSupportRingInnerRadius() + lugStemHalfLength; + if (verbosityLevel > 2) { + G4cout << __PRETTY_FUNCTION__ << " building support ring " << tag + << " positiveEnd=" << positiveEnd + << " ringHalfLength=" << ringHalfLength + << " inner/outer=" << tgt->stickmanSupportRingInnerRadius() + << "/" << tgt->stickmanSupportRingOuterRadius() + << " lugStemHalfLength=" << lugStemHalfLength + << " lugStemCenterRadius=" << lugStemCenterRadius + << G4endl; + } + + // ring body without lugs or center hole cutout + G4VSolid* supportRingSolid = reg.add(new G4Tubs( + "ProductionTargetSupportRingBase_" + tag, // + 0., // inner radius + tgt->stickmanSupportRingOuterRadius(), // outer radius + ringHalfLength, // half length + 0., // start angle + CLHEP::twopi)); // delta angle + + // cutout for the center of the ring + // such construction avoids potential issue that the lug solid extends within the + // inner radius of the ring, which can cause misrepresentation of the geometry + G4Tubs* ringCenterCutout = reg.add(new G4Tubs( + "ProductionTargetSupportRingCenterCutout_" + tag, // + 0., // inner radius + tgt->stickmanSupportRingInnerRadius(), // outer radius (avoid overlap with lug stem) + ringHalfLength + stickmanMagicOffset, // half length (add small offset to avoid precision issues) + 0., // start angle + CLHEP::twopi)); // delta angle + + // lugs at three fin angles. + G4Tubs* lugSolid = reg.add(new G4Tubs( + "ProductionTargetSupportRingLug_" + tag, // + 0., // inner radius + tgt->supportRingLugOuterRadius(), // outer radius + ringHalfLength, // half length + 0., // start angle + CLHEP::twopi)); // delta angle + + // lug stem + G4Box* lugStemSolid = reg.add(new G4Box( + "ProductionTargetSupportRingLugStem_" + tag, // + lugStemHalfLength, // half length in x + tgt->supportRingLugOuterRadius(), // half length in y + ringHalfLength)); // half length in z + + for (int ithFin = 0; ithFin < tgt->nStickmanFins(); ++ithFin) { + // the lugs are aligned with the fins, so get the fin angle for the current fin to calculate the position of the lugs + const double finAngle = tgt->plateFinAngle(ithFin); + const CLHEP::Hep3Vector lugStemShift( // rotate the lug stem shift by the fin angle to get the correct position in x and y + lugStemCenterRadius * std::cos(finAngle), + lugStemCenterRadius * std::sin(finAngle), + 0.); + const CLHEP::Hep3Vector lugShift( // rotate the lug shift by the fin angle + tgt->plateCenterToLugCenter() * std::cos(finAngle), + tgt->plateCenterToLugCenter() * std::sin(finAngle), + 0.); + if (verbosityLevel > 2) { + G4cout << __PRETTY_FUNCTION__ << " support ring " << tag + << " fin=" << ithFin + << " angle(deg)=" << finAngle/CLHEP::degree + << " lugStemShift=" << lugStemShift + << " lugShift=" << lugShift + << G4endl; + } + // first union the lug stem to the ring + supportRingSolid = reg.add(new G4UnionSolid( + "ProductionTargetSupportRingLugStemUnion_" + tag + "_" + std::to_string(ithFin), + supportRingSolid, // base solid + lugStemSolid, // union solid + stickmanFinRotations.at(ithFin), // passive rotation + lugStemShift)); // translation + // then union the lug to the ring+stem + supportRingSolid = reg.add(new G4UnionSolid( + "ProductionTargetSupportRingLugUnion_" + tag + "_" + std::to_string(ithFin), + supportRingSolid, // base solid + lugSolid, // union solid + nullptr, // passive rotation + lugShift)); // translation + + // add fillet to the lug if applicable + // The fillet is added by unioning with an extrusion + // of a triangle with vertices at (0,0), + // (sqrt((supportRingOuterRadius+filletRadius)^2-(supportRingLugOuterRadius+filletRadius)^2), supportRingLugOuterRadius+filletRadius), + // (sqrt((supportRingOuterRadius+filletRadius)^2-(supportRingLugOuterRadius+filletRadius)^2), -supportRingLugOuterRadius-filletRadius), + // extruded along z by stickmanSupportRingLength starting at the end of the target, with no twist and scale of 1.0 at both ends. + if (tgt->addFilletToSupportRingLug()) { + const double filletRadius = tgt->supportRingLugFilletRadius(); + const double filletX = std::sqrt( + std::pow(tgt->stickmanSupportRingOuterRadius() + filletRadius, 2) + - std::pow(tgt->supportRingLugOuterRadius() + filletRadius, 2)); + if (verbosityLevel > 2) { + G4cout << __PRETTY_FUNCTION__ << " support ring " << tag + << " fin=" << ithFin + << " filletRadius=" << filletRadius + << " filletX=" << filletX + << G4endl; } - VolumeInfo ringWithCutoutPositive(name,ringTranslation,prodTargetMotherInfo.centerInWorld); - ringWithCutoutPositive.solid = ringWithCutoutSolid.at(tgt->nHaymanFins() - 1); // what does this = mean? - finishNesting(ringWithCutoutPositive - ,prodTargetSupportRingMaterial - ,rotRing - ,ringTranslation - ,prodTargetMotherInfo.logical - ,finCopyNumber - ,G4Colour::White() - ,"ProductionTarget" - ,verbosityLevel>1); + // get the vertices for the extrusion + std::vector filletVertices = { + G4TwoVector(0., 0.), + G4TwoVector(filletX, tgt->supportRingLugOuterRadius() + filletRadius), + G4TwoVector(filletX, -tgt->supportRingLugOuterRadius() - filletRadius) + }; + + // create the extrusion solid for the fillet + G4ExtrudedSolid* filletExtrusion = reg.add(new G4ExtrudedSolid( + "ProductionTargetSupportRingLugFilletExtrusion_" + tag + "_" + std::to_string(ithFin), + filletVertices, // vertices + ringHalfLength, // half length in z + G4TwoVector(0., 0.), // shift of polygon center at -half length + 1.0, // scale of polygon at -half length + G4TwoVector(0., 0.), // shift of polygon center at +half length + 1.0)); // scale of polygon at +half length + // create the cutout cylinder for the fillet + G4Tubs* filletCutout = reg.add(new G4Tubs( + "ProductionTargetSupportRingLugFilletCutout_" + tag + "_" + std::to_string(ithFin), + 0., // inner radius + filletRadius + stickmanMagicOffset, // outer radius (add small offset) + ringHalfLength + stickmanMagicOffset, // half length in z (add small offset) + 0., // start angle + CLHEP::twopi)); // segment angle is full circle + // subtract the cutouts from the extrusion + G4VSolid* filletUpper = reg.add(new G4SubtractionSolid( + "ProductionTargetSupportRingLugFilletUpper_" + tag + "_" + std::to_string(ithFin), + filletExtrusion, // base solid + filletCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(filletX, + tgt->supportRingLugOuterRadius() + filletRadius, + 0.))); // translation of cutout solid + G4VSolid* filletSolid = reg.add(new G4SubtractionSolid( + "ProductionTargetSupportRingLugFilletLower_" + tag + "_" + std::to_string(ithFin), + filletUpper, // base solid + filletCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(filletX, + -tgt->supportRingLugOuterRadius() - filletRadius, + 0.))); // translation of cutout solid + + // add filletSolid to the union with the support ring + supportRingSolid = reg.add(new G4UnionSolid( + "ProductionTargetSupportRingFilletUnion_" + tag + "_" + std::to_string(ithFin), + supportRingSolid, // base solid + filletSolid, // union solid + stickmanFinRotations.at(ithFin), // passive rotation + CLHEP::Hep3Vector(0., 0., 0.))); // no translation + } // end of fillet construction for the lug + } // end of loop over fins to add lugs and fillets to the support ring + + // add cutout for the center hole + supportRingSolid = reg.add(new G4SubtractionSolid( + "ProductionTargetSupportRingLugFilletRingCutout_" + tag, + supportRingSolid, // base solid + ringCenterCutout, // cutout solid + nullptr, // rotation + CLHEP::Hep3Vector(0., 0., 0.))); // translation + + // add cutouts for the through holes if applicable + if (tgt->addCutoutToSupportRing()) { + const double cutoutTiltAngle = tgt->supportRingCutoutTilt(); + const double supportRingWallThickness = tgt->stickmanSupportRingOuterRadius() - tgt->stickmanSupportRingInnerRadius(); + // cutout center is placed at halfway between the inner and outer radius + const double cutoutCenterRadius = 0.5*(tgt->stickmanSupportRingInnerRadius() + tgt->stickmanSupportRingOuterRadius()); + // this is z position in the local ring frame + const double localCutoutZ = (positiveEnd ? 1. : -1.) * + (tgt->supportRingCutoutOffset() - ringHalfLength + supportRingWallThickness * 0.5 * std::tan(cutoutTiltAngle) ); + // cutout needs to go through the whole wall thickness + const double cutoutHalfLength = 0.5 * supportRingWallThickness / std::cos(cutoutTiltAngle) + + tgt->supportRingCutoutInnerRadius() * std::tan(cutoutTiltAngle) + + stickmanMagicOffset; // add small offset to avoid precision issues + if (verbosityLevel > 2) { + G4cout << __PRETTY_FUNCTION__ << " support ring " << tag + << " cutoutTilt(deg)=" << cutoutTiltAngle/CLHEP::degree + << " wallThickness=" << supportRingWallThickness + << " cutoutCenterRadius=" << cutoutCenterRadius + << " localCutoutZ=" << localCutoutZ + << " cutoutHalfLength=" << cutoutHalfLength + << G4endl; + } + // create the cutout solid for the through hole + G4Tubs* supportRingCutout = reg.add(new G4Tubs( + "ProductionTargetSupportRingCutout_" + tag, + 0., // inner radius + tgt->supportRingCutoutInnerRadius(), // outer radius + cutoutHalfLength, // half length in z + 0., // start angle + CLHEP::twopi)); // delta angle + const CLHEP::Hep3Vector localZAxis(0., 0., 1.); + + // loop over the holes and subtract the cutout solid + for (int ithCutout = 0; ithCutout < tgt->nSupportRingCutouts(); ++ithCutout) { + const double cutoutPhi = tgt->supportRingCutoutAngle(ithCutout); + const CLHEP::Hep3Vector radialAxis(std::cos(cutoutPhi), std::sin(cutoutPhi), 0.); // pointing from ring center radially outward + const CLHEP::Hep3Vector cutoutAxis = std::cos(cutoutTiltAngle) * radialAxis // vector pointing outward along the cutout axis + + (positiveEnd ? -1. : 1.) * std::sin(cutoutTiltAngle) * localZAxis; + G4RotationMatrix* cutoutRotation = reg.add(G4RotationMatrix( + CLHEP::HepRotation(cutoutAxis.cross(localZAxis), cutoutAxis.angle(localZAxis)))); // this rotates the cutout from the local z axis to the cutout axis passively + const CLHEP::Hep3Vector cutoutCenter( + cutoutCenterRadius * std::cos(cutoutPhi), + cutoutCenterRadius * std::sin(cutoutPhi), + localCutoutZ); + if (verbosityLevel > 2) { + G4cout << __PRETTY_FUNCTION__ << " support ring " << tag + << " cutout=" << ithCutout + << " phi(deg)=" << cutoutPhi/CLHEP::degree + << " center=" << cutoutCenter + << " axis=" << cutoutAxis + << G4endl; + } + // make this cutout solid + supportRingSolid = reg.add(new G4SubtractionSolid( + "ProductionTargetSupportRingCutoutSubtraction_" + tag + "_" + std::to_string(ithCutout), + supportRingSolid, // base solid + supportRingCutout, // cutout solid + cutoutRotation, // rotation to align the cutout with the cutout axis + cutoutCenter)); // translation + } // end of loop over cutouts for through holes + } // end of cutout construction for through holes + + return supportRingSolid; + }; // end of lambda to build the end support ring solid + + // now add the two end rings + const double supportRingCenterZ = tgt->halfStickmanLength() - tgt->stickmanSupportRingLength()/2.; + const CLHEP::Hep3Vector supportRingCenterNegZInTargetFrame(0., 0., -supportRingCenterZ); + const CLHEP::Hep3Vector supportRingCenterPosZInTargetFrame(0., 0., supportRingCenterZ); + const CLHEP::Hep3Vector supportRingCenterNegZ = + tgt->productionTargetRotation().inverse() * supportRingCenterNegZInTargetFrame; + const CLHEP::Hep3Vector supportRingCenterPosZ = + tgt->productionTargetRotation().inverse() * supportRingCenterPosZInTargetFrame; + if (verbosityLevel > 2) { + G4cout << __PRETTY_FUNCTION__ << " support ring centers target-frame neg/pos=" + << supportRingCenterNegZInTargetFrame << " / " + << supportRingCenterPosZInTargetFrame + << " mother-frame neg/pos=" + << supportRingCenterNegZ << " / " + << supportRingCenterPosZ + << G4endl; + } + // add the negative z end ring first + VolumeInfo supportRingNegZInfo("ProductionTargetSupportRing_Upstream", + supportRingCenterNegZ, + prodTargetMotherInfo.centerInWorld); + supportRingNegZInfo.solid = buildStickmanSupportRingSolid("Upstream", false); + finishNesting(supportRingNegZInfo, // volume info as defined above + prodTargetSupportRingMaterial, // material for the support ring + &tgt->productionTargetRotation(), // rotation + supportRingCenterNegZ, // translation + prodTargetMotherInfo.logical, // parent logical volume + 0, // copy number + G4Colour::Green(), + "ProductionTarget", + verbosityLevel>1); + // add the positive z end ring + VolumeInfo supportRingPosZInfo("ProductionTargetSupportRing_Downstream", + supportRingCenterPosZ, + prodTargetMotherInfo.centerInWorld); + supportRingPosZInfo.solid = buildStickmanSupportRingSolid("Downstream", true); + finishNesting(supportRingPosZInfo, // volume info as defined above + prodTargetSupportRingMaterial, // material for the support ring + &tgt->productionTargetRotation(), // rotation + supportRingCenterPosZ, // translation + prodTargetMotherInfo.logical, // parent logical volume + 1, // copy number + G4Colour::Green(), + "ProductionTarget", + verbosityLevel>1); + if (verbosityLevel > 0) { + G4cout << __PRETTY_FUNCTION__ << " end rings" + << " negZ center=" << supportRingCenterNegZ + << " posZ center=" << supportRingCenterPosZ + << " material=" << tgt->stickmanSupportRingMaterial() + << G4endl; + } - } //end if(ithSection == (numberOfSections -1) ) - } //end for(int ithSection...) - - //Add support structures for the production target - if(tgt->supportsBuild()) { - G4Material* suppWheelMaterial = findMaterialOrThrow(tgt->supportWheelMaterial()); - G4ThreeVector localWheelCenter(0.0,0.0,0.0); //no offset - double suppWheelParams[] = {tgt->supportWheelRIn(), tgt->supportWheelROut(), tgt->supportWheelHL()}; - //create the volume info for the support wheel+rods - VolumeInfo suppWheelInfo( "ProductionTargetSupportWheel", localWheelCenter, prodTargetMotherInfo.centerInMu2e()); - suppWheelInfo.solid = new G4Tubs("ProductionTargetSupportWheel_wheel", suppWheelParams[0], suppWheelParams[1], - suppWheelParams[2], 0., CLHEP::twopi); - // suppWheelParams, - // suppWheelMaterial, - // 0, - // localWheelCenter, - // prodTargetMotherInfo, - // 0, - // G4Colour::Gray(), - // "PS" - // ); - - // add spokes // - - //spoke info - const int nspokesperside = tgt->nSpokesPerSide(); - G4Material* spokeMaterial = findMaterialOrThrow(tgt->spokeMaterial()); - //target info - double rTarget = tgt->supportRingOuterRadius(); //radius of the support ring to attach to - double zTarget = tgt->halfHaymanLength(); //where along the target to attach - double smallGap = 0.001; //for adding small offsets to avoid overlaps due to precision - //initialize parameter vectors - //features on wheel - const vector supportWheelFeatureAngles = tgt->supportWheelFeatureAngles(); - const vector supportWheelFeatureArcs = tgt->supportWheelFeatureArcs (); - const vector supportWheelFeatureRIns = tgt->supportWheelFeatureRIns (); - //support rods in wheel - const vector supportWheelRodHL = tgt->supportWheelRodHL (); - const vector supportWheelRodOffset = tgt->supportWheelRodOffset (); - const vector supportWheelRodRadius = tgt->supportWheelRodRadius (); - const vector supportWheelRodRadialOffset = tgt->supportWheelRodRadialOffset(); - const vector supportWheelRodWireOffsetD = tgt->supportWheelRodWireOffsetD (); - const vector supportWheelRodWireOffsetU = tgt->supportWheelRodWireOffsetU (); - const vector supportWheelRodAngles = tgt->supportWheelRodAngles (); - //spoke (support wire) angles - const vector spokeTargetAnglesD = tgt->spokeTargetAnglesD(); - const vector spokeTargetAnglesU = tgt->spokeTargetAnglesU(); - if(verbosityLevel > 0) - std::cout << "Printing information about production target supports:\n"; - - const double targetAngle = tgt->rotHaymanY(); //assume target angle is only in the x-z plane for supports - CLHEP::HepRotation* rodRot = reg.add(CLHEP::HepRotation(CLHEP::HepRotation::IDENTITY)); - rodRot->rotateY(-1.*targetAngle); - - for(int istream = 0; istream < 2; ++istream) { - for(int ispoke = 0; ispoke < nspokesperside; ++ispoke) { - const double wheelAngle = supportWheelRodAngles[ispoke]*CLHEP::degree; - //get angle of the support rod on the wheel and the angle on the target the wire connects to - const double targetWireAngle = (istream == 0) ? spokeTargetAnglesD[ispoke]*CLHEP::degree - : spokeTargetAnglesU[ispoke]*CLHEP::degree; - double rWheel = supportWheelRodRadialOffset[ispoke]; // radius of the wheel to attach to - CLHEP::Hep3Vector rodCenter(rWheel*cos(wheelAngle), rWheel*sin(wheelAngle), 0.); - const double rodOffset = supportWheelRodOffset[ispoke]; - rodCenter += CLHEP::Hep3Vector(sin(targetAngle)*rodOffset, 0., cos(targetAngle)*rodOffset); - if(istream == 0) { //only do once - //add the features near the support rods in the bicycle wheel - const double featureAngle = supportWheelFeatureAngles[ispoke]*CLHEP::degree; //angle of feature center - const double featureArc = supportWheelFeatureArcs[ispoke]*CLHEP::degree; //width in angle - const double featureRIn = supportWheelFeatureRIns[ispoke]; //inner radius of feature - const double featureROut = tgt->supportWheelRIn() + smallGap; //ensure they overlap for union - // double featureR = (featureRIn + featureROut)/2.; //radius of feature center - // CLHEP::Hep3Vector featureCenter(featureR*cos(featureAngle), featureR*sin(featureAngle), 0.); - CLHEP::Hep3Vector featureCenter(localWheelCenter); //center is wheel center - double featureParams[] = {featureRIn, featureROut, tgt->supportWheelHL(), featureAngle - featureArc/2. /*phi0*/, featureArc /*dphi*/}; - G4Tubs* featureTubs = new G4Tubs("ProductionTargetSupportFeature_" +std::to_string(ispoke), - featureParams[0], featureParams[1], featureParams[2], featureParams[3], featureParams[4]); - suppWheelInfo.solid = new G4UnionSolid("ProductionTargetSupportWheelFeature_union_"+std::to_string(ispoke), - suppWheelInfo.solid, featureTubs, 0, featureCenter); - //add the support rod to the wheel - G4Tubs* rodTubs = new G4Tubs("ProductionTargetSupportRod_" +std::to_string(ispoke), - 0., supportWheelRodRadius[ispoke], supportWheelRodHL[ispoke], 0., CLHEP::twopi); - suppWheelInfo.solid = new G4UnionSolid("ProductionTargetSupportWheelRod_union_"+std::to_string(ispoke), - suppWheelInfo.solid, rodTubs, rodRot, rodCenter); - } - const int side = (1-2*istream); //+1 or -1 - //info about wire connection - //get end of the rod on this side - CLHEP::Hep3Vector rodAxis(sin(targetAngle), 0., cos(targetAngle)); - CLHEP::Hep3Vector wheelPos(rodCenter); - wheelPos += side*supportWheelRodHL[ispoke]*rodAxis; - //translate from rod center to edge - CLHEP::Hep3Vector rodCenterToWire(cos(wheelAngle)*cos(targetAngle), - sin(wheelAngle)*cos(targetAngle), - -cos(wheelAngle)*sin(targetAngle)); - wheelPos -= supportWheelRodRadius[ispoke]*rodCenterToWire; - double zWireRodOffset = (istream == 0) ? supportWheelRodWireOffsetD[ispoke] : supportWheelRodWireOffsetU[ispoke]; - wheelPos -= side*zWireRodOffset*rodAxis; - - //get wire position on target - CLHEP::Hep3Vector targetPos(rTarget*cos(targetWireAngle), rTarget*sin(targetWireAngle), side*zTarget); - targetPos = tgt->productionTargetRotation().inverse()*targetPos; //rotate from target frame to mother frame - CLHEP::Hep3Vector spokeAxis((wheelPos-targetPos).unit()); - CLHEP::Hep3Vector targetAxis(0.,0.,side); - targetAxis = tgt->productionTargetRotation().inverse()*targetAxis; - CLHEP::Hep3Vector zax(0.,0.,1.); - if(verbosityLevel > 0) - std::cout << "istream " << istream << " ispoke " << ispoke << std::endl - << "target pos " << targetPos << "\nwheel pos " << wheelPos << std::endl - << "Target axis " << targetAxis << "\nSpoke axis " << spokeAxis << std::endl - << "Rod axis " << rodAxis << "\nRod center to wire axis " << rodCenterToWire << std::endl; - //to remove overlaps where the wire connects, need angle of wire and surface connecting to - //remove overlap at target - double wireTargetAngle = targetAxis.angle(-1.*spokeAxis); - double deltaLength = (abs(tan(wireTargetAngle)) > 1.e-6) ? abs(tgt->spokeRadius()/tan(wireTargetAngle)) : 0.; //give up if ~paralle - targetPos += (deltaLength+0.1)*spokeAxis; //subtract off the length - if(verbosityLevel > 0) - std::cout << "wire target angle " << wireTargetAngle << " delta L " << deltaLength - << " target pos " << targetPos <spokeRadius()); - wheelPos -= (deltaLength+1.)*spokeAxis; - if(verbosityLevel > 0) - std::cout << "wire rod angle " << wireRodAngle << " delta L " << deltaLength - << " wheel pos " << wheelPos <spokeRadius(), 0.5*spokeLength); - CLHEP::HepRotation* spokeRot = reg.add(CLHEP::HepRotation(spokeAxis.cross(zax), spokeAxis.angle(zax))); - std::stringstream spokeName; - spokeName << "ProductionTargetSpokeWire_" ; - if(istream == 0) - spokeName << "Downstream_"; - else - spokeName << "Upstream_"; - spokeName << ispoke; - - VolumeInfo spokeInfo = nestTubs( spokeName.str(), - spokeParams, - spokeMaterial, - spokeRot, - spokeCenter, - prodTargetMotherInfo, - 0, - G4Colour::Gray(), - "PS" - ); - - } //end spokes loop - } //end stream loop - finishNesting(suppWheelInfo, - suppWheelMaterial, - 0, - localWheelCenter, - prodTargetMotherInfo.logical, - 0, - G4Colour::Gray(), - "PS" - ); + //Code following the construction of hayman_v_2_0 to add support structures. The Stickman support wheel + //is the same on design as the Hayman supports, but some errors and inaccuracies in the geometry implementations + //are fixed. + //Keeping it separate in the two branches instead of extracting it to a single helper function may seem redundant + //but is intentional. The Hayman geometry will exist as a benchmark in the future and in principle will not be + //changed or used. So it makes sense to keep it frozen and separated from the Stickman implementation, which may still + //see changes and iterations. Such "duplication" allows the two constructions to exist independently without affecting + //each other. + + //Add support structures for the production target + if(tgt->supportsBuild()) { + G4Material* suppWheelMaterial = findMaterialOrThrow(tgt->supportWheelMaterial()); + G4ThreeVector localWheelCenter(0.0,0.0,0.0); //no offset + double suppWheelParams[] = {tgt->supportWheelRIn(), tgt->supportWheelROut(), tgt->supportWheelHL()}; + //create the volume info for the support wheel+rods + VolumeInfo suppWheelInfo( "ProductionTargetSupportWheel", localWheelCenter, prodTargetMotherInfo.centerInMu2e()); + suppWheelInfo.solid = new G4Tubs("ProductionTargetSupportWheel_wheel", suppWheelParams[0], suppWheelParams[1], + suppWheelParams[2], 0., CLHEP::twopi); + // suppWheelParams, + // suppWheelMaterial, + // 0, + // localWheelCenter, + // prodTargetMotherInfo, + // 0, + // G4Colour::Gray(), + // "PS" + // ); + + // add spokes // + + //spoke info + const int nspokesperside = tgt->nSpokesPerSide(); + G4Material* spokeMaterial = findMaterialOrThrow(tgt->spokeMaterial()); + //target info + double rTarget = tgt->stickmanSupportRingOuterRadius(); //radius of the support ring to attach to + double zTarget = tgt->halfStickmanLength() + - tgt->stickmanSupportRingLength() + + tgt->supportRingCutoutOffset(); + //where along the target to attach, note the asymmetry + double smallGap = 0.001; //for adding small offsets to avoid overlaps due to precision + //initialize parameter vectors + //features on wheel + const vector supportWheelFeatureAngles = tgt->supportWheelFeatureAngles(); + const vector supportWheelFeatureArcs = tgt->supportWheelFeatureArcs (); + const vector supportWheelFeatureRIns = tgt->supportWheelFeatureRIns (); + //support rods in wheel + const vector supportWheelRodHL = tgt->supportWheelRodHL (); + const vector supportWheelRodOffset = tgt->supportWheelRodOffset (); + const vector supportWheelRodPinOffset = tgt->supportWheelRodPinOffset (); + const vector supportWheelRodRadius = tgt->supportWheelRodRadius (); + const vector supportWheelRodRadialOffset = tgt->supportWheelRodRadialOffset(); + const vector supportWheelRodWireOffsetD = tgt->supportWheelRodWireOffsetD (); + const vector supportWheelRodWireOffsetU = tgt->supportWheelRodWireOffsetU (); + const vector supportWheelRodAngles = tgt->supportWheelRodAngles (); + //spoke (support wire) angles + const vector spokeTargetAnglesD = tgt->spokeTargetAnglesD(); + const vector spokeTargetAnglesU = tgt->spokeTargetAnglesU(); + if(verbosityLevel > 0) + std::cout << "Printing information about production target supports:\n"; + + const double targetAngle = tgt->rotStickmanY(); //assume target angle is only in the x-z plane for supports + CLHEP::HepRotation* rodRot = reg.add(CLHEP::HepRotation(CLHEP::HepRotation::IDENTITY)); + rodRot->rotateY(-1.*targetAngle); + + for(int istream = 0; istream < 2; ++istream) { + for(int ispoke = 0; ispoke < nspokesperside; ++ispoke) { + const double wheelAngle = supportWheelRodAngles[ispoke]*CLHEP::degree; + //get angle of the support rod on the wheel and the angle on the target the wire connects to + const double targetWireAngle = (istream == 0) ? spokeTargetAnglesD[ispoke]*CLHEP::degree + : spokeTargetAnglesU[ispoke]*CLHEP::degree; + double rWheel = supportWheelRodRadialOffset[ispoke]; // radius of the wheel to attach to + CLHEP::Hep3Vector rodCenter(rWheel*std::cos(wheelAngle)/std::cos(targetAngle), + rWheel*std::sin(wheelAngle), + 0.); // rod plane projected to the wheel plane, non-orthogonal cut of a cylinder is an ellipse, + // so need to divide by cos(targetAngle) to get the correct position on the wheel plane + const double rodOffset = supportWheelRodOffset[ispoke] + + supportWheelRodPinOffset[ispoke]/std::cos(targetAngle); + rodCenter += CLHEP::Hep3Vector(std::sin(targetAngle)*rodOffset, 0., std::cos(targetAngle)*rodOffset); + if(istream == 0) { //only do once + //add the features near the support rods in the bicycle wheel + const double featureAngle = supportWheelFeatureAngles[ispoke]*CLHEP::degree; //angle of feature center + const double featureArc = supportWheelFeatureArcs[ispoke]*CLHEP::degree; //width in angle + const double featureRIn = supportWheelFeatureRIns[ispoke]; //inner radius of feature + const double featureROut = tgt->supportWheelRIn() + smallGap; //ensure they overlap for union + // double featureR = (featureRIn + featureROut)/2.; //radius of feature center + // CLHEP::Hep3Vector featureCenter(featureR*cos(featureAngle), featureR*sin(featureAngle), 0.); + CLHEP::Hep3Vector featureCenter(localWheelCenter); //center is wheel center + double featureParams[] = {featureRIn, featureROut, tgt->supportWheelHL(), featureAngle - featureArc/2. /*phi0*/, featureArc /*dphi*/}; + G4Tubs* featureTubs = new G4Tubs("ProductionTargetSupportFeature_" +std::to_string(ispoke), + featureParams[0], featureParams[1], featureParams[2], featureParams[3], featureParams[4]); + suppWheelInfo.solid = new G4UnionSolid("ProductionTargetSupportWheelFeature_union_"+std::to_string(ispoke), + suppWheelInfo.solid, featureTubs, 0, featureCenter); + //add the support rod to the wheel + G4Tubs* rodTubs = new G4Tubs("ProductionTargetSupportRod_" +std::to_string(ispoke), + 0., supportWheelRodRadius[ispoke], supportWheelRodHL[ispoke], 0., CLHEP::twopi); + suppWheelInfo.solid = new G4UnionSolid("ProductionTargetSupportWheelRod_union_"+std::to_string(ispoke), + suppWheelInfo.solid, rodTubs, rodRot, rodCenter); + } + const int side = (1-2*istream); //+1 or -1 + //info about wire connection + //get end of the rod on this side + CLHEP::Hep3Vector rodAxis(std::sin(targetAngle), 0., std::cos(targetAngle)); + CLHEP::Hep3Vector wheelPos(rodCenter); + wheelPos += side*supportWheelRodHL[ispoke]*rodAxis; + //translate from rod center to edge, this is vector from center to rod center rotated by target angle + CLHEP::Hep3Vector rodCenterToWire(std::cos(wheelAngle)*std::cos(targetAngle), + std::sin(wheelAngle), // rorated by target angle around y so y should be unchanged + -std::cos(wheelAngle)*std::sin(targetAngle)); + wheelPos -= supportWheelRodRadius[ispoke]*rodCenterToWire; + double zWireRodOffset = (istream == 0) ? supportWheelRodWireOffsetD[ispoke] : supportWheelRodWireOffsetU[ispoke]; + wheelPos -= side*zWireRodOffset*rodAxis; + + //get wire position on target + CLHEP::Hep3Vector targetPos(rTarget*std::cos(targetWireAngle), rTarget*std::sin(targetWireAngle), side*zTarget); + targetPos = tgt->productionTargetRotation().inverse()*targetPos; //rotate from target frame to mother frame + if (verbosityLevel > 1) + G4cout << __PRETTY_FUNCTION__ + << "Sanity check for spoke " << ispoke << " on stream " << istream << ":\n" + << "wheel pos " << wheelPos << "\ntarget pos " << targetPos << "\n" + << "spoke length " << (wheelPos-targetPos).mag() << G4endl; + CLHEP::Hep3Vector spokeAxis((wheelPos-targetPos).unit()); + CLHEP::Hep3Vector targetAxis(0.,0.,side); + targetAxis = tgt->productionTargetRotation().inverse()*targetAxis; + CLHEP::Hep3Vector zax(0.,0.,1.); + if(verbosityLevel > 0) + G4cout << __PRETTY_FUNCTION__ + << "istream " << istream << " ispoke " << ispoke << G4endl + << "target pos " << targetPos << "\nwheel pos " << wheelPos << G4endl + << "Target axis " << targetAxis << "\nSpoke axis " << spokeAxis << G4endl + << "Rod axis " << rodAxis << "\nRod center to wire axis " << rodCenterToWire << G4endl; + //to remove overlaps where the wire connects, need angle of wire and surface connecting to + //remove overlap at target + double wireTargetAngle = targetAxis.angle(-1.*spokeAxis); + double deltaLength = (std::abs(std::tan(wireTargetAngle)) > 1.e-6) ? std::abs(tgt->spokeRadius()/std::tan(wireTargetAngle)) : 0.; //give up if ~paralle + targetPos += (deltaLength+0.1)*spokeAxis; //subtract off the length + if(verbosityLevel > 0) + G4cout << __PRETTY_FUNCTION__ + << "wire target angle " << wireTargetAngle << " delta L " << deltaLength + << " target pos " << targetPos << G4endl; + + //next remove overlap at rod + // + double wireRodAngle = std::abs((rodCenterToWire).angle(spokeAxis)); + // deltaLength = std::abs(tan(wireRodAngle)/tgt->spokeRadius()); + // wheelPos -= (deltaLength+1.)*spokeAxis; + deltaLength = std::abs(tgt->spokeRadius()*std::cos(wireRodAngle)/std::sin(wireRodAngle)); // seems to be a typo in the hayman + wheelPos -= (deltaLength+0.1)*spokeAxis; + if(verbosityLevel > 0) + G4cout << __PRETTY_FUNCTION__ + << "wire rod angle " << wireRodAngle << " delta L " << deltaLength + << " wheel pos " << wheelPos << G4endl; + + CLHEP::Hep3Vector spokeCenter((wheelPos+targetPos)/2.); + double spokeLength = std::abs((wheelPos-targetPos).mag()); + TubsParams spokeParams(0., tgt->spokeRadius(), 0.5*spokeLength); + CLHEP::HepRotation* spokeRot = reg.add(CLHEP::HepRotation(spokeAxis.cross(zax), spokeAxis.angle(zax))); + std::stringstream spokeName; + spokeName << "ProductionTargetSpokeWire_" ; + if(istream == 0) + spokeName << "Downstream_"; + else + spokeName << "Upstream_"; + spokeName << ispoke; + + VolumeInfo spokeInfo = nestTubs( spokeName.str(), + spokeParams, + spokeMaterial, + spokeRot, + spokeCenter, + prodTargetMotherInfo, + 0, + G4Colour::Gray(), + "PS" + ); - } //end adding support structures - } //end ProductionTargetMaker::hayman_v_2_0 - } //end constructTargetPS + } //end spokes loop + } //end stream loop + finishNesting(suppWheelInfo, + suppWheelMaterial, + 0, + localWheelCenter, + prodTargetMotherInfo.logical, + 0, + G4Colour::Gray(), + "PS" + ); + + } //end adding support structures + } //end constructStickmanTarget } //end namespace mu2e - - // end Mu2eWorld::constructPS diff --git a/ProductionTargetGeom/inc/ProductionTarget.hh b/ProductionTargetGeom/inc/ProductionTarget.hh index 8ed4e3fac6..c94660563b 100644 --- a/ProductionTargetGeom/inc/ProductionTarget.hh +++ b/ProductionTargetGeom/inc/ProductionTarget.hh @@ -20,6 +20,100 @@ namespace mu2e { class ProductionTargetMaker; + // Parameter struct for transform and envelope parameters (Stickman) + struct StickmanEnvelopeParams { + double productionTargetMotherOuterRadius = 0.0; + double productionTargetMotherHalfLength = 0.0; + double rotStickmanX = 0.0; + double rotStickmanY = 0.0; + double rotStickmanZ = 0.0; + double halfStickmanLength = 0.0; + CLHEP::Hep3Vector stickmanProdTargetPosition; + std::string targetVacuumMaterial; + }; + + // Parameter struct for plate parameters (Stickman) + struct StickmanPlateParams { + int numberOfPlates = 0; + std::vector plateMaterial; + std::vector plateROut; + int nStickmanFins = 0; + std::vector plateFinAngles; + double plateFinOuterRadius = 0.0; + double plateFinWidth = 0.0; + double plateCenterToLugCenter = 0.0; + double plateLugInnerRadius = 0.0; + double plateLugOuterRadius = 0.0; + std::vector plateThickness; + std::vector plateLugThickness; + }; + + // Parameter struct for rod parameters (Stickman) + struct StickmanRodParams { + std::string rodMaterial; + double rodRadius = 0.0; + }; + + // Parameter struct for spacer parameters (Stickman) + struct StickmanSpacerParams { + std::string spacerMaterial; + double spacerHalfLength = 0.0; + double spacerOuterRadius = 0.0; + double spacerInnerRadius = 0.0; + }; + + // Parameter struct for support ring parameters (Stickman) + struct StickmanSupportRingParams { + std::string stickmanSupportRingMaterial; + double stickmanSupportRingLength = 0.0; + double stickmanSupportRingInnerRadius = 0.0; + double stickmanSupportRingOuterRadius = 0.0; + double supportRingLugOuterRadius = 0.0; + double supportRingCutoutOffset = 0.0; + }; + + // Parameter struct for Stickman configuration (additional parameters set after construction) + struct StickmanConfigParams { + // Plate fillet parameters + bool addFilletToPlateCore = false; + bool addFilletToPlateLug = false; + double plateFilletRadius = 0.0; + + // Support ring fillet parameters + bool addFilletToSupportRingLug = false; + double supportRingLugFilletRadius = 0.0; + + // Support ring cutout parameters + bool addCutoutToSupportRing = false; + int nSupportRingCutouts = 0; + std::vector supportRingCutoutAngles; + double supportRingCutoutInnerRadius = 0.0; + double supportRingCutoutTilt = 0.0; + + // Support wheel parameters + bool supportsBuild = false; + double supportWheelRIn = 0.0; + double supportWheelROut = 0.0; + double supportWheelHL = 0.0; + std::string supportWheelMaterial; + int nSpokesPerSide = 0; + std::vector supportWheelFeatureAngles; + std::vector supportWheelFeatureArcs; + std::vector supportWheelFeatureRIns; + std::vector supportWheelRodHL; + std::vector supportWheelRodOffset; + std::vector supportWheelRodPinOffset; + std::vector supportWheelRodRadius; + std::vector supportWheelRodRadialOffset; + std::vector supportWheelRodWireOffsetD; + std::vector supportWheelRodWireOffsetU; + std::vector supportWheelRodAngles; + std::vector spokeTargetAnglesD; + std::vector spokeTargetAnglesU; + double spokeRadius = 0.0; + std::string spokeMaterial; + }; + class ProductionTarget : virtual public Detector { public: @@ -123,6 +217,7 @@ namespace mu2e { //rods in the support wheel const std::vector& supportWheelRodHL () const {return _supportWheelRodHL ;} const std::vector& supportWheelRodOffset () const {return _supportWheelRodOffset ;} + const std::vector& supportWheelRodPinOffset () const {return _supportWheelRodPinOffset ;} //only used for Stickman_v_1_0 const std::vector& supportWheelRodRadius () const {return _supportWheelRodRadius ;} const std::vector& supportWheelRodRadialOffset() const {return _supportWheelRodRadialOffset;} const std::vector& supportWheelRodWireOffsetD () const {return _supportWheelRodWireOffsetD ;} @@ -137,27 +232,82 @@ namespace mu2e { double productionTargetMotherOuterRadius() const {return _productionTargetMotherOuterRadius;} double productionTargetMotherHalfLength() const {return _productionTargetMotherHalfLength;} + //accessors for Stickman_v_1_0 + std::string stickmanTargetType() const {return _stickmanTargetType;} + double halfStickmanLength() const {return _halfStickmanLength;} + int numberOfPlates() const {return _numberOfPlates;} + std::string plateMaterial(int i) const {return _plateMaterial.at(i);} + const std::vector& plateMaterial() const {return _plateMaterial;} + double plateROut(int i) const {return _plateROut.at(i);} + const std::vector& plateROut() const {return _plateROut;} + int nStickmanFins() const {return _nStickmanFins;} + double plateFinAngle(int i) const {return _plateFinAngles.at(i);} + const std::vector& plateFinAngles() const {return _plateFinAngles;} + double plateFinOuterRadius() const {return _plateFinOuterRadius;} + double plateFinWidth() const {return _plateFinWidth;} + double plateCenterToLugCenter() const {return _plateCenterToLugCenter;} + double plateLugInnerRadius() const {return _plateLugInnerRadius;} + double plateLugOuterRadius() const {return _plateLugOuterRadius;} + double plateThickness(int i) const {return _plateThickness.at(i);} + const std::vector& plateThickness() const {return _plateThickness;} + double plateLugThickness(int i) const {return _plateLugThickness.at(i);} + const std::vector& plateLugThickness() const {return _plateLugThickness;} + bool addFilletToPlateCore() const {return _addFilletToPlateCore;} + bool addFilletToPlateLug() const {return _addFilletToPlateLug;} + double plateFilletRadius() const {return _plateFilletRadius;} + std::string rodMaterial() const {return _rodMaterial;} + double rodRadius() const {return _rodRadius;} + double rodHalfLength() const {return _rodHalfLength;} + std::string spacerMaterial() const {return _spacerMaterial;} + double spacerHalfLength() const {return _spacerHalfLength;} + double spacerOuterRadius() const {return _spacerOuterRadius;} + double spacerInnerRadius() const {return _spacerInnerRadius;} + std::string stickmanSupportRingMaterial() const {return _stickmanSupportRingMaterial;} + double stickmanSupportRingLength() const {return _stickmanSupportRingLength;} + double stickmanSupportRingInnerRadius() const {return _stickmanSupportRingInnerRadius;} + double stickmanSupportRingOuterRadius() const {return _stickmanSupportRingOuterRadius;} + double supportRingLugOuterRadius() const {return _supportRingLugOuterRadius;} + bool addFilletToSupportRingLug() const {return _addFilletToSupportRingLug;} + double supportRingLugFilletRadius() const {return _supportRingLugFilletRadius;} + bool addCutoutToSupportRing() const {return _addCutoutToSupportRing;} + int nSupportRingCutouts() const {return _nSupportRingCutouts;} + double supportRingCutoutAngle(int i) const {return _supportRingCutoutAngles.at(i);} + const std::vector& supportRingCutoutAngles() const {return _supportRingCutoutAngles;} + double supportRingCutoutInnerRadius() const {return _supportRingCutoutInnerRadius;} + double supportRingCutoutTilt() const {return _supportRingCutoutTilt;} + double supportRingCutoutOffset() const {return _supportRingCutoutOffset;} + double rotStickmanX() const {return _rotStickmanX;} + double rotStickmanY() const {return _rotStickmanY;} + double rotStickmanZ() const {return _rotStickmanZ;} + CLHEP::Hep3Vector stickmanProdTargetPosition() const {return _stickmanProdTargetPosition;} + std::string hayman_v_2_0 = "Hayman_v_2_0"; std::string tier1 = "MDC2018"; + std::string stickman_v_1_0 = "Stickman_v_1_0"; CLHEP::Hep3Vector targetPositionByVersion() const { if (_haymanTargetType == hayman_v_2_0){ return _haymanProdTargetPosition;} - else if (_tier1TargetType == "MDC2018"){ + else if (_tier1TargetType == tier1){ return _prodTargetPosition;} + else if (_stickmanTargetType == stickman_v_1_0){ + return _stickmanProdTargetPosition;} else throw cet::exception("BADCONFIG") << "in ProductionTarget.hh, no valid target specified"<< std::endl; } double targetHalfLengthByVersion() const { if (_haymanTargetType == hayman_v_2_0){ return _halfHaymanLength;} - else if (_tier1TargetType == "MDC2018"){ + else if (_tier1TargetType == tier1){ return _halfLength;} + else if (_stickmanTargetType == stickman_v_1_0){ + return _halfStickmanLength;} else throw cet::exception("BADCONFIG") << "in ProductionTarget.hh, no valid target specified"<< std::endl; } - + // Configuration method for additional Stickman parameters + void configureStickman(const StickmanConfigParams& configParams); //---------------------------------------------------------------- @@ -205,6 +355,16 @@ namespace mu2e { ,double supportRingCutoutLength ); + ProductionTarget( + std::string stickmanTargetType + ,int version + ,const StickmanEnvelopeParams& envelopeParams + ,const StickmanPlateParams& plateParams + ,const StickmanRodParams& rodParams + ,const StickmanSpacerParams& spacerParams + ,const StickmanSupportRingParams& supportRingParams + ); + CLHEP::HepRotation _protonBeamRotation; // can't return by const ref if invert on the fly so need to store redundant data @@ -289,6 +449,7 @@ namespace mu2e { //parameters for rods in the support wheel std::vector _supportWheelRodHL ; //includes length through the wheel std::vector _supportWheelRodOffset ; //z offset with respect to the wheel + std::vector _supportWheelRodPinOffset ; //pinhole offset from support wheel center plane, used for Stickman only std::vector _supportWheelRodRadius ; //radius of the rod std::vector _supportWheelRodRadialOffset; //radius from the wheel center the rod is centered at std::vector _supportWheelRodWireOffsetD ; //z offset from the end of the rod the wire connects (downstream) @@ -302,6 +463,49 @@ namespace mu2e { std::vector _spokeTargetAnglesU; //angle about the target the wire connects to (upstream) double _spokeRadius ; //radius of the wire + // Stickman_v_1_0 parameters + std::string _stickmanTargetType; + double _halfStickmanLength; + double _rotStickmanX; + double _rotStickmanY; + double _rotStickmanZ; + CLHEP::Hep3Vector _stickmanProdTargetPosition; + int _numberOfPlates; + std::vector _plateMaterial; + std::vector _plateROut; + int _nStickmanFins; + std::vector _plateFinAngles; + double _plateFinOuterRadius; + double _plateFinWidth; + double _plateCenterToLugCenter; + double _plateLugInnerRadius; + double _plateLugOuterRadius; + std::vector _plateThickness; + std::vector _plateLugThickness; + bool _addFilletToPlateCore; + bool _addFilletToPlateLug; + double _plateFilletRadius; + std::string _rodMaterial; + double _rodRadius; + double _rodHalfLength; + std::string _spacerMaterial; + double _spacerHalfLength; + double _spacerOuterRadius; + double _spacerInnerRadius; + std::string _stickmanSupportRingMaterial; + double _stickmanSupportRingLength; + double _stickmanSupportRingInnerRadius; + double _stickmanSupportRingOuterRadius; + double _supportRingLugOuterRadius; + bool _addFilletToSupportRingLug; + double _supportRingLugFilletRadius; + bool _addCutoutToSupportRing; + int _nSupportRingCutouts; + std::vector _supportRingCutoutAngles; + double _supportRingCutoutInnerRadius; + double _supportRingCutoutTilt; + double _supportRingCutoutOffset; + // Needed for persistency template friend class art::Wrapper; ProductionTarget():_pHubsRgtParams(NULL), _pHubsLftParams(NULL) {} diff --git a/ProductionTargetGeom/src/ProductionTarget.cc b/ProductionTargetGeom/src/ProductionTarget.cc index 042fad66b4..4011b92eea 100644 --- a/ProductionTargetGeom/src/ProductionTarget.cc +++ b/ProductionTargetGeom/src/ProductionTarget.cc @@ -1,4 +1,6 @@ #include "Offline/ProductionTargetGeom/inc/ProductionTarget.hh" +#include "cetlib_except/exception.h" +#include namespace mu2e { ProductionTarget::ProductionTarget(std::string tier1TargetType, int version, double rOut, @@ -99,5 +101,181 @@ namespace mu2e { _prodTargetPosition = _haymanProdTargetPosition; } + ProductionTarget::ProductionTarget( + std::string stickmanTargetType + ,int version + ,const StickmanEnvelopeParams& envelopeParams + ,const StickmanPlateParams& plateParams + ,const StickmanRodParams& rodParams + ,const StickmanSpacerParams& spacerParams + ,const StickmanSupportRingParams& supportRingParams + ) + : _protonBeamRotation(CLHEP::HepRotation::IDENTITY) + ,_version(version) + ,_productionTargetMotherOuterRadius(envelopeParams.productionTargetMotherOuterRadius) + ,_productionTargetMotherHalfLength(envelopeParams.productionTargetMotherHalfLength) + ,_targetVacuumMaterial(envelopeParams.targetVacuumMaterial) + ,_stickmanTargetType(stickmanTargetType) + ,_halfStickmanLength(envelopeParams.halfStickmanLength) + ,_rotStickmanX(envelopeParams.rotStickmanX) + ,_rotStickmanY(envelopeParams.rotStickmanY) + ,_rotStickmanZ(envelopeParams.rotStickmanZ) + ,_stickmanProdTargetPosition(envelopeParams.stickmanProdTargetPosition) + ,_numberOfPlates(plateParams.numberOfPlates) + ,_plateMaterial(plateParams.plateMaterial) + ,_plateROut(plateParams.plateROut) + ,_nStickmanFins(plateParams.nStickmanFins) + ,_plateFinAngles(plateParams.plateFinAngles) + ,_plateFinOuterRadius(plateParams.plateFinOuterRadius) + ,_plateFinWidth(plateParams.plateFinWidth) + ,_plateCenterToLugCenter(plateParams.plateCenterToLugCenter) + ,_plateLugInnerRadius(plateParams.plateLugInnerRadius) + ,_plateLugOuterRadius(plateParams.plateLugOuterRadius) + ,_plateThickness(plateParams.plateThickness) + ,_plateLugThickness(plateParams.plateLugThickness) + ,_rodMaterial(rodParams.rodMaterial) + ,_rodRadius(rodParams.rodRadius) + ,_spacerMaterial(spacerParams.spacerMaterial) + ,_spacerHalfLength(spacerParams.spacerHalfLength) + ,_spacerOuterRadius(spacerParams.spacerOuterRadius) + ,_spacerInnerRadius(spacerParams.spacerInnerRadius) + ,_stickmanSupportRingMaterial(supportRingParams.stickmanSupportRingMaterial) + ,_stickmanSupportRingLength(supportRingParams.stickmanSupportRingLength) + ,_stickmanSupportRingInnerRadius(supportRingParams.stickmanSupportRingInnerRadius) + ,_stickmanSupportRingOuterRadius(supportRingParams.stickmanSupportRingOuterRadius) + ,_supportRingLugOuterRadius(supportRingParams.supportRingLugOuterRadius) + ,_supportRingCutoutOffset(supportRingParams.supportRingCutoutOffset) + { + if (_stickmanTargetType.empty()) { + throw cet::exception("GEOM") << "ProductionTarget: missing stickman target type"; + } + if (_productionTargetMotherOuterRadius <= 0.0) { + throw cet::exception("GEOM") << "ProductionTarget: invalid productionTargetMotherOuterRadius"; + } + if (_productionTargetMotherHalfLength <= 0.0) { + throw cet::exception("GEOM") << "ProductionTarget: invalid productionTargetMotherHalfLength"; + } + if (_halfStickmanLength <= 0.0) { + throw cet::exception("GEOM") << "ProductionTarget: invalid halfStickmanLength"; + } + if (_numberOfPlates <= 0) { + throw cet::exception("GEOM") << "ProductionTarget: numberOfPlates must be positive"; + } + if (_plateMaterial.size() != static_cast(_numberOfPlates)) { + throw cet::exception("GEOM") + << "ProductionTarget: targetPS_plateMaterial size mismatch: expected " + << _numberOfPlates << ", got " << _plateMaterial.size(); + } + if (_plateROut.size() != static_cast(_numberOfPlates)) { + throw cet::exception("GEOM") + << "ProductionTarget: targetPS_rOut size mismatch: expected " + << _numberOfPlates << ", got " << _plateROut.size(); + } + if (_plateThickness.size() != static_cast(_numberOfPlates)) { + throw cet::exception("GEOM") + << "ProductionTarget: targetPS_plateThickness size mismatch: expected " + << _numberOfPlates << ", got " << _plateThickness.size(); + } + if (_plateLugThickness.size() != static_cast(_numberOfPlates)) { + throw cet::exception("GEOM") + << "ProductionTarget: targetPS_plateLugThickness size mismatch: expected " + << _numberOfPlates << ", got " << _plateLugThickness.size(); + } + if (_nStickmanFins <= 0) { + throw cet::exception("GEOM") << "ProductionTarget: nStickmanFins must be positive"; + } + if (_plateFinAngles.size() != static_cast(_nStickmanFins)) { + throw cet::exception("GEOM") + << "ProductionTarget: targetPS_plateFinAngles size mismatch: expected " + << _nStickmanFins << ", got " << _plateFinAngles.size(); + } + if (_rodMaterial.empty()) { + throw cet::exception("GEOM") << "ProductionTarget: missing targetPS_rodMaterial"; + } + if (_rodRadius <= 0.0) { + throw cet::exception("GEOM") << "ProductionTarget: invalid targetPS_rodRadius"; + } + if (_spacerMaterial.empty()) { + throw cet::exception("GEOM") << "ProductionTarget: missing targetPS_spacerMaterial"; + } + if (_spacerHalfLength <= 0.0) { + throw cet::exception("GEOM") << "ProductionTarget: invalid targetPS_spacerHalfLength"; + } + if (_spacerOuterRadius <= 0.0) { + throw cet::exception("GEOM") << "ProductionTarget: invalid targetPS_spacerOuterRadius"; + } + if (_spacerInnerRadius < 0.0 || _spacerInnerRadius >= _spacerOuterRadius) { + throw cet::exception("GEOM") << "ProductionTarget: invalid spacer radii"; + } + if (_stickmanSupportRingMaterial.empty()) { + throw cet::exception("GEOM") << "ProductionTarget: missing targetPS_supportRingMaterial"; + } + if (_stickmanSupportRingLength <= 0.0) { + throw cet::exception("GEOM") << "ProductionTarget: invalid targetPS_supportRingLength"; + } + if (_stickmanSupportRingInnerRadius <= 0.0) { + throw cet::exception("GEOM") << "ProductionTarget: invalid targetPS_supportRingInnerRadius"; + } + if (_stickmanSupportRingOuterRadius <= 0.0) { + throw cet::exception("GEOM") << "ProductionTarget: invalid targetPS_supportRingOuterRadius"; + } + if (_supportRingLugOuterRadius <= 0.0) { + throw cet::exception("GEOM") << "ProductionTarget: invalid targetPS_supportRingLugOuterRadius"; + } + if (_supportRingCutoutOffset <= 0.0) { + throw cet::exception("GEOM") << "ProductionTarget: invalid targetPS_supportRingCutoutOffset"; + } + + // rod half length, actual rod is longer than this since it inserts into the end rings, but for the geometry reconstruction, that part will be treated as part of the end ring. + _rodHalfLength = std::accumulate(_plateLugThickness.begin(), _plateLugThickness.end(), 0.0) / 2.0 + 2 * _spacerHalfLength; + + // rotations + _protonBeamRotation.rotateX(envelopeParams.rotStickmanX).rotateY(envelopeParams.rotStickmanY).rotateZ(envelopeParams.rotStickmanZ); + _protonBeamInverseRotation = _protonBeamRotation.inverse(); // passive rotation matrix + _halfLength = _productionTargetMotherHalfLength; + _prodTargetPosition = _stickmanProdTargetPosition; + } + + void ProductionTarget::configureStickman(const StickmanConfigParams& configParams) { + // Configure plate fillet parameters + _addFilletToPlateCore = configParams.addFilletToPlateCore; + _addFilletToPlateLug = configParams.addFilletToPlateLug; + _plateFilletRadius = configParams.plateFilletRadius; + + // Configure support ring fillet parameters + _addFilletToSupportRingLug = configParams.addFilletToSupportRingLug; + _supportRingLugFilletRadius = configParams.supportRingLugFilletRadius; + + // Configure support ring cutout parameters + _addCutoutToSupportRing = configParams.addCutoutToSupportRing; + _nSupportRingCutouts = configParams.nSupportRingCutouts; + _supportRingCutoutAngles = configParams.supportRingCutoutAngles; + _supportRingCutoutInnerRadius = configParams.supportRingCutoutInnerRadius; + _supportRingCutoutTilt = configParams.supportRingCutoutTilt; + + // Configure support wheel parameters + _supportsBuild = configParams.supportsBuild; + _supportWheelRIn = configParams.supportWheelRIn; + _supportWheelROut = configParams.supportWheelROut; + _supportWheelHL = configParams.supportWheelHL; + _supportWheelMaterial = configParams.supportWheelMaterial; + _nSpokesPerSide = configParams.nSpokesPerSide; + _supportWheelFeatureAngles = configParams.supportWheelFeatureAngles; + _supportWheelFeatureArcs = configParams.supportWheelFeatureArcs; + _supportWheelFeatureRIns = configParams.supportWheelFeatureRIns; + _supportWheelRodHL = configParams.supportWheelRodHL; + _supportWheelRodOffset = configParams.supportWheelRodOffset; + _supportWheelRodPinOffset = configParams.supportWheelRodPinOffset; + _supportWheelRodRadius = configParams.supportWheelRodRadius; + _supportWheelRodRadialOffset = configParams.supportWheelRodRadialOffset; + _supportWheelRodWireOffsetD = configParams.supportWheelRodWireOffsetD; + _supportWheelRodWireOffsetU = configParams.supportWheelRodWireOffsetU; + _supportWheelRodAngles = configParams.supportWheelRodAngles; + _spokeTargetAnglesD = configParams.spokeTargetAnglesD; + _spokeTargetAnglesU = configParams.spokeTargetAnglesU; + _spokeRadius = configParams.spokeRadius; + _spokeMaterial = configParams.spokeMaterial; + } + } diff --git a/ProductionTargetGeom/src/SConscript b/ProductionTargetGeom/src/SConscript index 30c881054e..b0c342fb6e 100644 --- a/ProductionTargetGeom/src/SConscript +++ b/ProductionTargetGeom/src/SConscript @@ -13,7 +13,8 @@ helper=mu2e_helper(env); mainlib = helper.make_mainlib ( [ 'mu2e_GeomPrimitives', 'CLHEP', - 'Core' + 'Core', + 'cetlib_except', ] ) # This tells emacs to view this file in python mode.