diff --git a/DataProducts/inc/SurfaceId.hh b/DataProducts/inc/SurfaceId.hh index ae8fdddbd9..0c298133b3 100644 --- a/DataProducts/inc/SurfaceId.hh +++ b/DataProducts/inc/SurfaceId.hh @@ -18,6 +18,7 @@ namespace mu2e { TT_Front=0, TT_Mid, TT_Back, TT_Inner, TT_Outer, // tracker VD equivalents DS_Front=80, DS_Back, DS_Inner, DS_Outer, IPA_Legacy, //introduced for backwards compatibility with MDS1 + DS_CryoInner, DS_CryoOuter, DS_ShieldInner, DS_ShieldOuter, DS_Coil, IPA=90, IPA_Front, IPA_Back, OPA=95, TSDA, // Absorbers in the DS ST_Front=100,ST_Back, ST_Inner, ST_Outer, ST_Foils, ST_Wires, // stopping target bounding surfaces and components @@ -32,11 +33,13 @@ namespace mu2e { CRV_U =240, // CRV-Upstream CRV_D1=250, CRV_D2, CRV_D3, CRV_D4, // CRV-Downstream CRV_C1=260, CRV_C2, // CRV-Cryo-Outer - CRV_M1=270, CRV_M2, CRV_M3, CRV_M4, CRV_M5, CRV_M6, CRV_M7, CRV_M8 //CRV-Muon-Taggers + CRV_M1=270, CRV_M2, CRV_M3, CRV_M4, CRV_M5, CRV_M6, CRV_M7, CRV_M8, // CRV-Muon-Taggers (Mu2e/Offline PR #1864) + CRV_StrongBack=280, // CRV module Al strongback (tracker-side support plate) + DS_HatchConcrete=300 // detector-area hatch concrete block approximation }; // Update this counter whenever you add/remove surface IDs from the enum above. - static constexpr std::size_t nSurfaceIds = 55; + static constexpr std::size_t nSurfaceIds = 62; static std::string const& typeName(); static std::map const& names(); diff --git a/DataProducts/src/SurfaceId.cc b/DataProducts/src/SurfaceId.cc index 10de06ec48..e984a523bf 100644 --- a/DataProducts/src/SurfaceId.cc +++ b/DataProducts/src/SurfaceId.cc @@ -28,6 +28,11 @@ namespace mu2e { std::make_pair(SurfaceIdEnum::DS_Inner, "DS_Inner"), std::make_pair(SurfaceIdEnum::DS_Outer, "DS_Outer"), std::make_pair(SurfaceIdEnum::IPA_Legacy, "IPA_Legacy"), + std::make_pair(SurfaceIdEnum::DS_CryoInner, "DS_CryoInner"), + std::make_pair(SurfaceIdEnum::DS_CryoOuter, "DS_CryoOuter"), + std::make_pair(SurfaceIdEnum::DS_ShieldInner, "DS_ShieldInner"), + std::make_pair(SurfaceIdEnum::DS_ShieldOuter, "DS_ShieldOuter"), + std::make_pair(SurfaceIdEnum::DS_Coil, "DS_Coil"), std::make_pair(SurfaceIdEnum::IPA, "IPA"), std::make_pair(SurfaceIdEnum::IPA_Front, "IPA_Front"), std::make_pair(SurfaceIdEnum::IPA_Back, "IPA_Back"), @@ -71,7 +76,9 @@ namespace mu2e { std::make_pair(SurfaceIdEnum::CRV_M5, "CRV_M5"), std::make_pair(SurfaceIdEnum::CRV_M6, "CRV_M6"), std::make_pair(SurfaceIdEnum::CRV_M7, "CRV_M7"), - std::make_pair(SurfaceIdEnum::CRV_M8, "CRV_M8") + std::make_pair(SurfaceIdEnum::CRV_M8, "CRV_M8"), + std::make_pair(SurfaceIdEnum::CRV_StrongBack, "CRV_StrongBack"), + std::make_pair(SurfaceIdEnum::DS_HatchConcrete, "DS_HatchConcrete") }; constexpr std::size_t nSurfaceIdNames = sizeof(surfaceIdNames)/sizeof(surfaceIdNames[0]); diff --git a/ExternalShieldingGeom/inc/ExtShieldDownstream.hh b/ExternalShieldingGeom/inc/ExtShieldDownstream.hh index ed03a75b55..1ffb9ed47c 100644 --- a/ExternalShieldingGeom/inc/ExtShieldDownstream.hh +++ b/ExternalShieldingGeom/inc/ExtShieldDownstream.hh @@ -70,6 +70,13 @@ namespace mu2e { { return _notchLocations; } const std::vector >& getNotchDimensions() const { return _notchDimensions; } + // 1-based box-type number per block (the "type" from the SimpleConfig description). + const std::vector& getBoxTypes() const + { return _boxTypes; } + // Declared number of box types (ExtShieldDownstream.numberOfBoxTypes); may exceed the + // number of types that actually have blocks (some types can have zero blocks). + int getNumberOfBoxTypes() const + { return _numberOfBoxTypes; } private: @@ -91,7 +98,9 @@ namespace mu2e { const std::vector& oHole, const std::vector& iNotch, const std::vector& locNotch, - const std::vector >& locDims) + const std::vector >& locDims, + const std::vector& boxTypes, + int numberOfBoxTypes) : _extShieldOutlines(outlines), _extShieldLengths (lengths), _extShieldBoxTols (tols), @@ -107,7 +116,9 @@ namespace mu2e { _holeOrientations (oHole), _notchIndices (iNotch), _notchLocations (locNotch), - _notchDimensions (locDims) + _notchDimensions (locDims), + _boxTypes (boxTypes), + _numberOfBoxTypes (numberOfBoxTypes) { } // Or read back from persistent storage @@ -137,6 +148,9 @@ namespace mu2e { std::vector< int > _notchIndices; std::vector< CLHEP::Hep3Vector > _notchLocations; std::vector< std::vector< double > > _notchDimensions; + // 1-based box-type number per block (parallel to _centerPositions etc.) + std::vector< int > _boxTypes; + int _numberOfBoxTypes = 0; }; std::ostream& operator<<(std::ostream& os, const ExtShieldDownstream& upstr); diff --git a/GeometryService/inc/KinKalGeomMaker.hh b/GeometryService/inc/KinKalGeomMaker.hh index 4b9a4f8106..af8fd18549 100644 --- a/GeometryService/inc/KinKalGeomMaker.hh +++ b/GeometryService/inc/KinKalGeomMaker.hh @@ -6,16 +6,26 @@ // #include "Offline/KinKalGeom/inc/KinKalGeom.hh" #include "Offline/KinKalGeom/inc/KKMaterial.hh" +#include "CLHEP/Vector/ThreeVector.h" +#include namespace mu2e { + class DetectorSystem; class KinKalGeomMaker { public: - KinKalGeomMaker(int debug) : debug_(debug) {} + explicit KinKalGeomMaker(int debug) : debug_(debug) {} std::unique_ptr& makeKKG(); private: void makeTracker(); void makeDS(); void makeTarget(); void makeCRV(); + void makePassiveMaterials(); + // append one averaged rectangular concrete passive-material plane (see .cc) + void addConcretePlane(DetectorSystem const& det, int normalAxis, + CLHEP::Hep3Vector const& centerMu2e, double hw1, double hw2, + double halfThickness, std::string const& material); + // run-1 building hatch concrete when ExtShieldDownstream concrete is zeroed + void addBuildingHatchConcrete(DetectorSystem const& det); std::unique_ptr kkg_; int debug_ = 0; }; diff --git a/GeometryService/src/ExtShieldDownstreamMaker.cc b/GeometryService/src/ExtShieldDownstreamMaker.cc index 045e29bca9..86aa6db854 100644 --- a/GeometryService/src/ExtShieldDownstreamMaker.cc +++ b/GeometryService/src/ExtShieldDownstreamMaker.cc @@ -150,6 +150,7 @@ namespace mu2e { std::vector holeOri; std::vector notchLoc; std::vector > notchDim; + std::vector boxTypes; outlines .reserve(nBoxesTot); lengths .reserve(nBoxesTot); tols .reserve(nBoxesTot); @@ -160,6 +161,7 @@ namespace mu2e { nNotches .reserve(nBoxesTot); holeID .reserve(nBoxesTot); notchID .reserve(nBoxesTot); + boxTypes .reserve(nBoxesTot); // Helper variables for building lists of holes and notches int holeIdx = 0; @@ -172,6 +174,7 @@ namespace mu2e { lengths.push_back(lengthOfType[it]); // from the type-specified tols.push_back(tolsOfType[it]); // info filled above mats.push_back(materialOfType[it]); + boxTypes.push_back(it+1); // 1-based type number for this box // Location of the center of the box in Mu2e coords @@ -301,7 +304,9 @@ namespace mu2e { holeOri, notchID, notchLoc, - notchDim) + notchDim, + boxTypes, + nType) ); //---------------------------------------------------------------- diff --git a/GeometryService/src/KinKalGeomMaker.cc b/GeometryService/src/KinKalGeomMaker.cc index a4af3b909f..37191bb987 100644 --- a/GeometryService/src/KinKalGeomMaker.cc +++ b/GeometryService/src/KinKalGeomMaker.cc @@ -3,8 +3,13 @@ // Original author: Dave Brown (LBNL) 4/2026 // #include "Offline/GeometryService/inc/KinKalGeomMaker.hh" +#include "Offline/GeometryService/inc/GeometryService.hh" #include "Offline/TrackerGeom/inc/Tracker.hh" #include "Offline/CosmicRayShieldGeom/inc/CosmicRayShield.hh" +#include "Offline/CosmicRayShieldGeom/inc/CRSScintillatorShield.hh" +#include "Offline/CosmicRayShieldGeom/inc/CRSScintillatorModule.hh" +#include "Offline/ExternalShieldingGeom/inc/ExtShieldDownstream.hh" +#include "Offline/Mu2eHallGeom/inc/Mu2eHall.hh" #include "Offline/GeometryService/inc/GeomHandle.hh" #include "Offline/KinKalGeom/inc/KinKalGeom.hh" #include "Offline/KinKalGeom/inc/Tracker.hh" @@ -14,9 +19,15 @@ #include "Offline/BeamlineGeom/inc/Beamline.hh" #include "Offline/GeometryService/inc/DetectorSystem.hh" #include "Offline/DetectorSolenoidGeom/inc/DetectorSolenoid.hh" +#include "CLHEP/Vector/ThreeVector.h" +#include "CLHEP/Vector/TwoVector.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" #include "cetlib_except/exception.h" #include #include +#include +#include +#include namespace mu2e { using KinKal::VEC3; @@ -42,6 +53,7 @@ namespace mu2e { makeDS(); makeTarget(); makeCRV(); + makePassiveMaterials(); return kkg_; } @@ -123,7 +135,71 @@ namespace mu2e { auto outer= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-1482),ds->rOut2(),ds->halfLength()); // bounding surfaces auto front= std::make_shared(outer->frontDisk()); auto back= std::make_shared(outer->backDisk()); + KKGeom::DetectorSolenoid::MaterialCylinderCollection materialCylinders; + auto toDetectorZ = [&det](double zmu2e) { + return VEC3(det->toDetector(CLHEP::Hep3Vector(0.0,0.0,zmu2e))).Z(); + }; + auto addMaterialCylinder = [this,&materialCylinders,&toDetectorZ](SurfaceId const& sid, + double rin, double rout, double zcenter, double halfLength, std::string const& material) { + auto cylinder = std::make_shared(VEC3(0.0,0.0,1.0), + VEC3(0.0,0.0,toDetectorZ(zcenter)),0.5*(rin+rout),halfLength); + materialCylinders.emplace_back(sid,cylinder,material,rout-rin); + kkg_->map_.emplace(std::make_pair(sid,std::static_pointer_cast(cylinder))); + }; + // KinKal DetMaterial name for the DS thermal shield (aluminium). The G4 geometry calls this + // material "G4_Al"; here we use a DS-specific MatEnv material so the passive plane self-documents. + std::string const dsShieldMaterial = "DSShieldAl"; + addMaterialCylinder(SurfaceId(SurfaceIdEnum::DS_CryoInner),ds->rIn1(),ds->rIn2(),ds->position().z(),ds->halfLength(),ds->material()); + addMaterialCylinder(SurfaceId(SurfaceIdEnum::DS_CryoOuter),ds->rOut1(),ds->rOut2(),ds->position().z(),ds->halfLength(),ds->material()); + addMaterialCylinder(SurfaceId(SurfaceIdEnum::DS_ShieldInner),ds->shield_rIn1(),ds->shield_rIn2(),ds->position().z(),ds->shield_halfLength(),dsShieldMaterial); + addMaterialCylinder(SurfaceId(SurfaceIdEnum::DS_ShieldOuter),ds->shield_rOut1(),ds->shield_rOut2(),ds->position().z(),ds->shield_halfLength(),dsShieldMaterial); + // Coils: collapse the 11 short coil cylinders into ONE long cylinder spanning the coil z-range. + // NB: the DSCoilMix density (TrackerConditions/data/MaterialsList.data) assumes the effective thickness below; + // both follow the run-2 coil geometry. + { + size_t const ncoils = static_cast(ds->nCoils()); + double zmin = 1e9, zmax = -1e9, tsum = 0.0, rsum = 0.0, zlentot = 0.0; + for(size_t icoil = 0; icoil < ncoils; icoil++){ + double const zlen = ds->coil_zLength().at(icoil); + double const z0 = ds->coil_zPosition().at(icoil); + double const rout = ds->coil_rOut().at(icoil); + zmin = std::min(zmin, z0); zmax = std::max(zmax, z0 + zlen); + tsum += zlen * (rout - ds->coil_rIn()); // z-weighted radial thickness + rsum += zlen * 0.5*(ds->coil_rIn() + rout); // z-weighted mid radius + zlentot += zlen; + } + double const teff = tsum / zlentot; // effective radial thickness (~28.08 mm) + double const reff = rsum / zlentot; // effective mid radius + addMaterialCylinder(SurfaceId(SurfaceIdEnum::DS_Coil), + reff - 0.5*teff, reff + 0.5*teff, 0.5*(zmin + zmax), 0.5*(zmax - zmin), "DSCoilMix"); + } + // DS cryostat + thermal-shield END WALLS (z-normal annuli). Tracks bound for the upstream/ + // downstream CRV sectors (CRV-U / CRV-D) leave the solenoid through its z-faces and cross the + // cryostat and shield end walls; the radial shells above only catch tracks that exit radially. + // Geometry from DetectorSolenoid: (shield_)endWallHalfLength(). + { + double const dsz = ds->position().z(); + auto addEndWall = [this,&toDetectorZ](SurfaceId const& sid, double zmu2e, double rin, + double rout, double halfThick, std::string const& material) { + auto annulus = std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0), + VEC3(0.0,0.0,toDetectorZ(zmu2e)),rin,rout); + kkg_->passiveMaterialPlanes_.emplace_back(sid,annulus,material,halfThick); + kkg_->map_.emplace(std::make_pair(sid,std::static_pointer_cast(annulus))); + }; + double const cOff = ds->halfLength() - ds->endWallHalfLength(); // cryo wall center inset + double const sOff = ds->shield_halfLength() - ds->shield_endWallHalfLength(); // shield wall center inset + // cryostat end walls (StainlessSteel) + addEndWall(SurfaceId(SurfaceIdEnum::DS_Front,1), dsz-cOff, ds->rIn1(), ds->rOut2(), + ds->endWallHalfLength(), ds->material()); + addEndWall(SurfaceId(SurfaceIdEnum::DS_Back,1), dsz+cOff, ds->rIn1(), ds->rOut2(), + ds->endWallHalfLength(), ds->material()); + // thermal-shield end walls (aluminium; see dsShieldMaterial above) + addEndWall(SurfaceId(SurfaceIdEnum::DS_Front,2), dsz-sOff, ds->shield_rIn1(), ds->shield_rOut2(), + ds->shield_endWallHalfLength(), dsShieldMaterial); + addEndWall(SurfaceId(SurfaceIdEnum::DS_Back,2), dsz+sOff, ds->shield_rIn1(), ds->shield_rOut2(), + ds->shield_endWallHalfLength(), dsShieldMaterial); + } // hard-coded for now auto ipa= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-2770),300.0,500.0); @@ -142,7 +218,7 @@ namespace mu2e { kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::OPA),std::static_pointer_cast(opa))); kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::TSDA),std::static_pointer_cast(tsda))); - kkg_->ds_ = std::make_unique(inner, outer, front, back, ipa, ipafront, ipaback, opa, tsda); + kkg_->ds_ = std::make_unique(inner, outer, front, back, ipa, ipafront, ipaback, opa, tsda, materialCylinders); } void KinKalGeomMaker::makeTarget() { @@ -177,6 +253,16 @@ namespace mu2e { GeomHandle det; auto const& shields = CRS->getCRSScintillatorShields(); std::vector sectors; + // Strongback (aluminium support plate) thickness from the CRV geometry + double strongBackThickness = 0.0; + for(auto const& shield : shields){ + int const td = shield.getFirstBar().getBarDetail().getThicknessDirection(); + for(int im = 0; im < shield.nModules(); ++im){ + auto const& mod = shield.getModule(im); + if(mod.nAluminumSheets() > 0){ strongBackThickness = 2.0*mod.getAluminumSheet(0).getHalfLengths().at(td); break; } + } + if(strongBackThickness != 0.0) break; + } // loop over the shields (= sectors) for (auto const& shield : shields) { // @@ -261,6 +347,17 @@ namespace mu2e { sector.sname_ = shield.getName(); sector.sector_ = std::make_shared(wdir,udir,midpoint,uhw,vhw); sector.whw_ = whw; + // Each CRV module carries a ~12.7 mm aluminium strongback (support plate) on the tracker-facing side + // of the scintillator stack. Thickness (strongBackThickness, above) and material now come from the + // CRV geometry rather than crs.strongBackThickness / crs.aluminumSheetMaterialName. + // KinKal DetMaterial name for the strongback (aluminium). The CRV geometry's aluminium-sheet + // material is the G4 name "G4_Al"; we use a CRV-specific MatEnv material so the plane self-documents. + std::string const strongBackMaterial = "CRVStrongbackAl"; + auto const sbcenter = midpoint - (whw + 0.5*strongBackThickness)*wdir; + auto const sbplane = std::make_shared(wdir,udir,sbcenter,uhw,vhw); + SurfaceId sbsid(SurfaceIdEnum::CRV_StrongBack, static_cast(kkg_->passiveMaterialPlanes_.size())); + kkg_->passiveMaterialPlanes_.emplace_back(sbsid, sbplane, strongBackMaterial, 0.5*strongBackThickness); + kkg_->map_.emplace(std::make_pair(sbsid, std::static_pointer_cast(sbplane))); sectors.push_back(sector); } // sort the sectors according to their transverse distance (largest first), to optimize searching for downward going tracks. @@ -282,4 +379,368 @@ namespace mu2e { isect++; } } + + // Append one averaged rectangular concrete passive-material plane for an ExtShield region. + // 'halfThickness' is the half material depth a track crosses along the normal + // (used for the energy-loss/scattering path length). All planes use the broad + // DS_HatchConcrete SurfaceId; a running index distinguishes the individual planes. + void KinKalGeomMaker::addConcretePlane(DetectorSystem const& det, + int normalAxis, CLHEP::Hep3Vector const& centerMu2e, + double hw1, double hw2, double halfThickness, std::string const& material) { + // Build the plane's local frame. The Rectangle ctor takes (normal,udir,center, + // uHalfLen,vHalfLen) with vdir = normal x udir (right-handed). We pick udir/vdir + // as the two global axes orthogonal to the normal so the half-extents map cleanly. + VEC3 norm, udir; + double uhw, vhw; + switch(normalAxis) { + case 0: // x-normal (side walls): u=y, v=z + norm = VEC3(1.0,0.0,0.0); udir = VEC3(0.0,1.0,0.0); uhw = hw1; vhw = hw2; break; + case 1: // y-normal (roof): u=x, v=z (matches the Type-2 roof convention below) + norm = VEC3(0.0,1.0,0.0); udir = VEC3(1.0,0.0,0.0); uhw = hw1; vhw = hw2; break; + case 2: // z-normal (downstream endcap): u=x, v=y + norm = VEC3(0.0,0.0,1.0); udir = VEC3(1.0,0.0,0.0); uhw = hw1; vhw = hw2; break; + default: + throw cet::exception("Service") << "invalid concrete-plane normal axis " << normalAxis << std::endl; + } + // Split the crossing into two Gauss-point planes along the normal, at center +/- halfThickness/sqrt3, + // each carrying material half-depth 'halfThickness' (so the two together = the full slab depth + // 2*halfThickness). This integrates the bulk multiple scattering with a 2-point Gauss rule (exact + // for a uniform slab: reproduces the continuous angle-position covariance) while preserving the total + // crossed material, so energy loss is correct. + double const invSqrt3 = 1.0/std::sqrt(3.0); + auto const center0 = VEC3(det.toDetector(centerMu2e)); + for(int iplane = 0; iplane < 2; ++iplane) { + double const sign = iplane == 0 ? -1.0 : 1.0; + auto const center = center0 + (sign*halfThickness*invSqrt3)*norm; + auto const plane = std::make_shared(norm, udir, center, uhw, vhw); + SurfaceId sid(SurfaceIdEnum::DS_HatchConcrete, + static_cast(kkg_->passiveMaterialPlanes_.size())); + kkg_->passiveMaterialPlanes_.emplace_back(sid, plane, material, halfThickness); + kkg_->map_.emplace(std::make_pair(sid, std::static_pointer_cast(plane))); + } + } + + void KinKalGeomMaker::makePassiveMaterials() { + GeomHandle det; + // External shielding (concrete) geometry now comes from the ExtShieldDownstream GeometryService object + art::ServiceHandle geomService; + if(!geomService->hasElement()) return; + GeomHandle esd; + if(esd->getNumberOfBoxTypes() <= 0) return; + // Group block indices by their 1-based type number. The type identity is what SimpleConfig + // supplied and what outline/material/orientation alone cannot recover (several types share + // identical geometry); getBoxTypes() restores it. Vectors below are parallel per block. + auto const& esdTypes = esd->getBoxTypes(); + auto const& esdOutlines = esd->getOutlines(); // per block: list of {u,v} vertices (mm) + auto const& esdLengths = esd->getLengths(); // per block: HALF length along W (mm) + auto const& esdMats = esd->getMaterialNames(); + auto const& esdCenters = esd->getCentersOfBoxes(); + std::map> boxesByType; + for(std::size_t i = 0; i < esdTypes.size(); ++i) boxesByType[esdTypes[i]].push_back(i); + auto firstBox = [&](int t)->long { + auto it = boxesByType.find(t); + return (it == boxesByType.end() || it->second.empty()) ? -1L : static_cast(it->second.front()); + }; + auto nBoxOf = [&](int t){ auto it = boxesByType.find(t); return it == boxesByType.end() ? 0 : static_cast(it->second.size()); }; + auto materialOf = [&](int t){ long const b = firstBox(t); return b >= 0 ? esdMats[b] : std::string(); }; + // Mu2e-frame center {x,y,z} of the ibox-th (1-based) block of a type (mirrors centerType{t}Box{ibox}). + auto centerOf = [&](int t, int ibox)->std::vector { + auto const& c = esdCenters[boxesByType[t][ibox-1]]; + return { c.x(), c.y(), c.z() }; + }; + // --- ExtShieldDownstream box dimensions. A box "type" is an extruded U-V outline of half-length + // esdLengths; the helpers return its half-extents. The {U,V,W=length} -> {x,y,z} orientation per + // region is not carried here, so it stays explicit at the call sites (as before). + // Helpers return 0 for absent types so an optional region is simply skipped (its nBox guard fires). + auto vertsOf = [&](int t, bool uaxis){ + std::vector v; long const b = firstBox(t); + if(b >= 0) for(auto const& uv : esdOutlines[b]) v.push_back(uaxis ? uv[0] : uv[1]); + return v; + }; + auto maxAbs = [](std::vector const& v){ double m=0.0; for(double x : v) m=std::max(m,std::fabs(x)); return m; }; + auto spanHalf= [](std::vector const& v){ if(v.empty()) return 0.0; auto e=std::minmax_element(v.begin(),v.end()); return 0.5*(*e.second - *e.first); }; + auto uHalf = [&](int t){ return maxAbs(vertsOf(t,true)); }; // outline U half-width (max |U|) + auto vHalf = [&](int t){ return maxAbs(vertsOf(t,false)); }; // outline V half-height (max |V|) + auto vSpanH = [&](int t){ return spanHalf(vertsOf(t,false)); }; // V half-extent (handles 0-based outlines) + auto lenHalf = [&](int t){ long const b = firstBox(t); return b >= 0 ? esdLengths[b] : 0.0; }; // already half + // V-weighted mean half-thickness of a T cross-section along its U (thin) axis: used where a + // T-block is collapsed to a single slab whose normal is U (the side walls), giving the + // energy-loss-correct mean depth without per-component Gauss planes. Reproduces the run-2 449.7. + auto tMeanUHalf = [&](int t){ + auto uV=vertsOf(t,true), vV=vertsOf(t,false); + if(uV.empty() || vV.empty()) return 0.0; + double const uo=maxAbs(uV); double ui=uo; + for(double u : uV) if(u > 0.0 && u < uo-1.0) ui=std::min(ui,u); // inner (stem) U half-width + auto const ve=std::minmax_element(vV.begin(),vV.end()); + double const vmin=*ve.first, vmax=*ve.second; double vsh=0.5*(vmin+vmax); + for(double v : vV) if(v > vmin+1.0 && v < vmax-1.0){ vsh=v; break; } // shoulder V + double const cb=vsh-vmin, st=vmax-vsh, span=vmax-vmin; // crossbar / stem V extents + return span>0.0 ? (uo*cb + ui*st)/span : uo; + }; + + // Model the ExtShieldDownstream Type 2 regular concrete roof T-blocks. + // These blocks sit between the DS outer cryostat (y~1328 mm) and the CRV top sectors T1/T2 (y~2653 mm). + // With orientations "010"/"012" (both map U->z, V->y, W->x): + // U (±uhw_outer mm) -> ±z footprint of the wide crossbar per block + // V range [vmin, vmax] -> y; vshoulder separates crossbar from stem + // W (length/2 = xhw mm) -> x (transverse, same for both components) + // + // The T cross-section is decomposed into two non-overlapping rectangles: + // Crossbar: V in [vmin, vshoulder], U in [-uhw_outer, +uhw_outer] + // Stem : V in [vshoulder, vmax], U in [-uhw_inner, +uhw_inner] + // Each is modelled with two Gauss-point planes so tracks through the crossbar- + // only region (the ~2/3 of z-width outside the stem) correctly see half the + // concrete compared to tracks through the full T height. + // Wrapped in an IIFE so a degenerate Type-2 (no boxes / empty outline) skips only this block + // via `return`, without abandoning the side walls / roof / endcap / hatch handled below. + [&]{ + int const nType2 = nBoxOf(2); + if(nType2 <= 0) return; + double const invSqrt3 = 1.0/std::sqrt(3.0); + std::vector uVerts = vertsOf(2,true), vVerts = vertsOf(2,false); + if(uVerts.empty() || vVerts.empty()) return; + + auto const uExt = std::minmax_element(uVerts.begin(), uVerts.end()); + auto const vExt = std::minmax_element(vVerts.begin(), vVerts.end()); + double const xhw = lenHalf(2); + auto const material = materialOf(2); + + double const vmin_t = *vExt.first; // = -452.2 mm (bottom of crossbar) + double const vmax_t = *vExt.second; // = +452.2 mm (top of stem) + double const uhw_outer = *uExt.second; // = 680.8 mm (crossbar z half-width per block) + + // Inner U half-width (stem): smallest positive U vertex that differs from uhw_outer + double uhw_inner = uhw_outer; + for(double u : uVerts) + if(u > 0 && u < uhw_outer - 1.0) uhw_inner = std::min(uhw_inner, u); + + // Shoulder V (where T width steps inward): V value strictly between vmin and vmax + double vshoulder = 0.0; + for(double v : vVerts) + if(v > vmin_t + 1.0 && v < vmax_t - 1.0) { vshoulder = v; break; } + + // Crossbar component: V in [vmin_t, vshoulder] + double const yhalf_cb = 0.5*(vshoulder - vmin_t); // = 223.6 mm + double const vcenter_cb = 0.5*(vmin_t + vshoulder); // = -228.6 mm (local V offset) + + // Stem component: V in [vshoulder, vmax_t] + double const yhalf_stem = 0.5*(vmax_t - vshoulder); // = 228.6 mm + double const vcenter_stem = 0.5*(vshoulder + vmax_t); // = +223.6 mm (local V offset) + + // Accumulate x/y center (common to both components) and z bounding boxes + // (different per component because stem is narrower in z than crossbar). + double xcenter = 0.0, ycenter = 0.0; + double zmin_cb = 1.0e9, zmax_cb = -1.0e9; + double zmin_st = 1.0e9, zmax_st = -1.0e9; + int nfound = 0; + for(int ibox = 1; ibox <= nType2; ++ibox) { + std::vector ctr = centerOf(2, ibox); + xcenter += ctr[0]; + ycenter += ctr[1]; + zmin_cb = std::min(zmin_cb, ctr[2] - uhw_outer); + zmax_cb = std::max(zmax_cb, ctr[2] + uhw_outer); + zmin_st = std::min(zmin_st, ctr[2] - uhw_inner); + zmax_st = std::max(zmax_st, ctr[2] + uhw_inner); + ++nfound; + } + if(nfound == 0) return; + xcenter /= nfound; + ycenter /= nfound; + + // Global Mu2e Y centers for each T component (ycenter = mean box center = 2006.1 mm) + double const ycenter_cb = ycenter + vcenter_cb; // ≈ 1777.5 mm + double const ycenter_stem = ycenter + vcenter_stem; // ≈ 2229.7 mm + + double const zcenters[2] = { 0.5*(zmin_cb + zmax_cb), 0.5*(zmin_st + zmax_st) }; + double const zhws[2] = { 0.5*(zmax_cb - zmin_cb), 0.5*(zmax_st - zmin_st) }; + double const yctrs[2] = { ycenter_cb, ycenter_stem }; + double const yhalfs[2] = { yhalf_cb, yhalf_stem }; + + // Four planes total: two Gauss-point planes for each T component. + for(int icomp = 0; icomp < 2; ++icomp) { + for(int iplane = 0; iplane < 2; ++iplane) { + double const ysign = iplane == 0 ? -1.0 : 1.0; + auto const center = VEC3(det->toDetector( + CLHEP::Hep3Vector(xcenter, yctrs[icomp] + ysign*yhalfs[icomp]*invSqrt3, zcenters[icomp]))); + auto const plane = std::make_shared( + VEC3(0.0,1.0,0.0), VEC3(1.0,0.0,0.0), center, xhw, zhws[icomp]); + SurfaceId sid(SurfaceIdEnum::DS_HatchConcrete, + static_cast(kkg_->passiveMaterialPlanes_.size())); + kkg_->passiveMaterialPlanes_.emplace_back(sid, plane, material, yhalfs[icomp]); + kkg_->map_.emplace(std::make_pair(sid, std::static_pointer_cast(plane))); + } + } + }(); + + // Collapse all type- ExtShieldDownstream concrete boxes of a region into one + // averaged plane. 'normalAxis' is the slab normal (0=x,1=y,2=z); the in-plane + // half-extents are (center span)/2 + the block half-size along each in-plane + // axis; 'halfThickness' is supplied by the caller (block depth along normal). + // Boxes whose center lies far from the region cluster (different sub-region, + // e.g. the type-1 north vs south wall) are separated by an explicit center cut. + auto buildAveragedRegion = [&](int boxType, int normalAxis, + double halfU_block, double halfV_block, double halfThickness, + std::function const&)> select) { + int const nbox = nBoxOf(boxType); + if(nbox <= 0) return; + auto const material = materialOf(boxType); + // axis indices orthogonal to the normal, in increasing order + int a1 = (normalAxis == 0) ? 1 : 0; + int a2 = (normalAxis == 2) ? 1 : 2; + double nsum = 0.0; // accumulate the normal-axis center (slab position) + double min1 = 1e9, max1 = -1e9; // in-plane axis 1 center range + double min2 = 1e9, max2 = -1e9; // in-plane axis 2 center range + int nfound = 0; + for(int ibox = 1; ibox <= nbox; ++ibox) { + std::vector ctr = centerOf(boxType, ibox); + if(select && !select(ctr)) continue; + nsum += ctr[normalAxis]; + min1 = std::min(min1, ctr[a1]); max1 = std::max(max1, ctr[a1]); + min2 = std::min(min2, ctr[a2]); max2 = std::max(max2, ctr[a2]); + ++nfound; + } + if(nfound == 0) return; + // assemble the center from the per-axis values (normal axis = mean slab + // position; in-plane axes = midpoint of the center bounding box) + double cc[3] = {0.0,0.0,0.0}; + cc[normalAxis] = nsum / nfound; + cc[a1] = 0.5*(min1 + max1); + cc[a2] = 0.5*(min2 + max2); + CLHEP::Hep3Vector center(cc[0],cc[1],cc[2]); + double const hw1 = 0.5*(max1 - min1) + halfU_block; + double const hw2 = 0.5*(max2 - min2) + halfV_block; + this->addConcretePlane(*det, normalAxis, center, hw1, hw2, halfThickness, material); + }; + + // --- Side walls. Type 1 regular-concrete wall T-blocks. Geometry from + // ExtShieldDownstream (outlineType1 U/V verts + lengthType1): the T cross-section + // lies in the U-V plane, extruded over lengthType1; for the walls the orientation + // maps U->x (wall thickness), V->z (footprint), length->y (vertical extent). + // U verts +/-680.8 / +/-223.6 -> crossbar 1361.6 mm thick, stem 447.2 mm thick + // V verts +/-452.2 (span 904.4); crossbar occupies V in [-452.2,-5] = 447.2 mm, + // stem occupies V in [-5, 452.2] = 457.2 mm + // The wall is one x-normal slab per side whose x half-thickness is the z-weighted + // mean of the T thickness, so the mean crossed concrete (energy loss) is exact + // without four Gauss planes per wall: + // 2*xHalf = (1361.6*447.2 + 447.2*457.2)/904.4 = 899.3 mm -> xHalfThick = 449.7 + // Box centers split into the north wall (x ~ -1897) and south wall (x ~ -5911) + // about the DS axis (x = -3904); each is emitted as a single slab. + { + double const yHalfBlock = lenHalf(1); // length (W) -> y vertical half-height + double const zHalfBlock = vSpanH(1); // V outline -> z footprint (452.2 in run-2) + double const xHalfThick = tMeanUHalf(1); // U-outline T-mean -> x half-thickness (449.7 in run-2) + // Separate the two symmetric side walls at the DS axis = the mean x of the type-1 box + // centers (-3904 in run-2), derived from config so it tracks the geometry instead of a literal. + double splitX = 0.0; int nWall = 0; + int const nBox1 = nBoxOf(1); + for(int ibox = 1; ibox <= nBox1; ++ibox) { + std::vector c = centerOf(1, ibox); + splitX += c[0]; ++nWall; + } + if(nWall > 0) { + splitX /= nWall; + buildAveragedRegion(1, 0, yHalfBlock, zHalfBlock, xHalfThick, + [splitX](std::vector const& c){ return c[0] > splitX; }); // north + buildAveragedRegion(1, 0, yHalfBlock, zHalfBlock, xHalfThick, + [splitX](std::vector const& c){ return c[0] < splitX; }); // south + } + } + + // --- Other roof concrete types (same y-normal plane family as Type-2). These + // sit at y ~ 2006 between the DS and the top CRV. We add + // them as averaged y-normal slabs using the roof T half-height (the V outline + // ~452.2 mm maps to the vertical thickness for the roof). Types: 4 (barite roof + // T, 6 boxes), 9 (barite roof L, 2), 10 (concrete roof filler, 1), 23 (barite + // roof T around ST, 2), 29 (endcap top-back, 1), 30 (endcap top-front, 1). + // For these the wide U crossbar (680.8) and the W length (4918 for the long + // roof types) span x and z; we use the conservative full-width block half-sizes + // so the averaged footprint covers the tiled blocks. y half-thickness = 452.2. + { + double const roofYhalfThick = vHalf(2); // V outline half -> roof vertical depth (452.2) + double const roofXhalfBlock = lenHalf(2); // length/2 spans x (~2459) + double const roofZhalfBlock = uHalf(2); // U crossbar half spans z (680.8) + // These slabs borrow Type-2's dimensions; if Type-2 is zeroed but these types aren't, skip + // rather than emit zero-thickness planes (a partial config the run-1 nBoxOf(2)<=0 path misses). + if(roofYhalfThick > 0.0) { + // type 4: barite roof T blocks (upstream) + buildAveragedRegion(4, 1, roofXhalfBlock, roofZhalfBlock, roofYhalfThick, nullptr); + // type 23: barite roof T blocks around the stopping target + buildAveragedRegion(23, 1, roofXhalfBlock, roofZhalfBlock, roofYhalfThick, nullptr); + // type 9: barite roof L blocks on the ends of the upstream roof + buildAveragedRegion(9, 1, roofXhalfBlock, roofZhalfBlock, roofYhalfThick, nullptr); + // type 10: concrete roof filler block + buildAveragedRegion(10, 1, roofXhalfBlock, roofZhalfBlock, roofYhalfThick, nullptr); + // types 29/30: new-endcap top back / front roof pieces. y-normal slabs; their own wide U + // outline spans x and their short length spans z, modelled at the roof slab depth. + buildAveragedRegion(29, 1, uHalf(29), lenHalf(29), roofYhalfThick, nullptr); + buildAveragedRegion(30, 1, uHalf(30), lenHalf(30), roofYhalfThick, nullptr); + } + } + + // --- Downstream endcap concrete (z-normal). Types 27 (3 boxes) and 28 (1 box) + // form the wall closing the DS<->CRV gap at z ~ 18542, beyond the DS, at the + // CRV-D sectors. They have large x/y outlines (U up to 2463.8, V up to 3246.2, + // the V outline runs 0..3246.2 so the per-block V half-extent is ~1623 mm) and + // a short z length (~1117 / ~907 mm), i.e. a slab normal to z. We emit one + // averaged z-normal plane per type, with x half from the U outline, y half from + // the centers' bounding box plus the block V half, and z half-thickness from + // length/2. NOTE: the asymmetric (0-based) V outline means the precise y + // placement of the face would need the per-box "300" rotation to be decoded; + // here we cover the whole stacked x-y face conservatively, which is adequate for + // the broad energy-loss correction this region needs for CRV-D-bound tracks. + { + // Endcap face (z-normal): the type-27 U outline -> x, its V outline span -> y; the slab + // z half-thickness is each type's own length/2. (Type 28 shares the type-27 x-y face, as before.) + double const ecXhalf = uHalf(27); // U outline -> x half (2463.8) + double const ecYhalf = vSpanH(27); // V outline span -> y half (~1623) + // Both endcap types use Type-27's x-y face; skip if Type-27 is zeroed so Type-28 doesn't + // emit a zero-area plane. + if(ecXhalf > 0.0 && ecYhalf > 0.0) { + // type 27: 3 stacked endcap blocks (cover most of the x-y face) + buildAveragedRegion(27, 2, ecXhalf, ecYhalf, lenHalf(27), nullptr); + // type 28: endcap bottom block + buildAveragedRegion(28, 2, ecXhalf, ecYhalf, lenHalf(28), nullptr); + } + } + + // Run-1 geometries zero the ExtShieldDownstream concrete T-blocks but retain the building + // dsAreaHatchBlock concrete (CONCRETE_MARS) between tracker and CRV-T. Extracted geometries + // have no ExtShieldDownstream object and returned above. Reaching here means the shield is + // declared (numberOfBoxTypes>0); the run-1 fallback fires when its roof T-blocks are zeroed. + if(nBoxOf(2) <= 0) { + addBuildingHatchConcrete(*det); + } + } + + void KinKalGeomMaker::addBuildingHatchConcrete(DetectorSystem const& det) { + // The dsAreaHatchBlock concrete is a building ExtrudedSolid served by Mu2eHall, replacing the + // former building.dsArea.hatchblock.* / yOfFloorSurface SimpleConfig reads. The solid's offset + // already folds the floor-surface Y into its Mu2e-origin offset. The 2D outline lies in the + // (x, z) plane: vertex x -> Mu2e x, vertex y -> Mu2e z; the slab normal is y (a roof plane). + GeomHandle hall; + auto const& solids = hall->getBldgSolids(); + auto spanHalf = [](std::vector const& v, bool xaxis) { + double mn = 1e30, mx = -1e30; + for(auto const& p : v){ double const c = xaxis ? p.x() : p.y(); mn = std::min(mn,c); mx = std::max(mx,c); } + return 0.5*(mx - mn); + }; + auto mid = [](std::vector const& v, bool xaxis) { + double mn = 1e30, mx = -1e30; + for(auto const& p : v){ double const c = xaxis ? p.x() : p.y(); mn = std::min(mn,c); mx = std::max(mx,c); } + return 0.5*(mn + mx); + }; + auto addSlab = [&](std::string const& name) { + auto const it = solids.find(name); + if(it == solids.end()) return; + auto const& solid = it->second; + auto const& verts = solid.getVertices(); + if(verts.empty()) return; + auto const& off = solid.getOffsetFromMu2eOrigin(); + CLHEP::Hep3Vector center(off.x() + mid(verts,true), off.y(), off.z() + mid(verts,false)); + addConcretePlane(det, 1, center, spanHalf(verts,true), spanHalf(verts,false), + solid.getYhalfThickness(), solid.getMaterial()); + }; + addSlab("dsAreaHatchBlock"); + addSlab("dsAreaHatchBlockE"); + } } diff --git a/KinKalGeom/inc/DetectorSolenoid.hh b/KinKalGeom/inc/DetectorSolenoid.hh index 5fd98a55db..c697d3b27e 100644 --- a/KinKalGeom/inc/DetectorSolenoid.hh +++ b/KinKalGeom/inc/DetectorSolenoid.hh @@ -8,7 +8,9 @@ #include "KinKal/Geometry/Frustrum.hh" #include "KinKal/Geometry/Annulus.hh" #include "KinKal/Geometry/Disk.hh" +#include "Offline/DataProducts/inc/SurfaceId.hh" #include +#include #include namespace mu2e { namespace KKGeom { @@ -18,12 +20,22 @@ namespace mu2e { using FruPtr = std::shared_ptr; using DiskPtr = std::shared_ptr; using AnnPtr = std::shared_ptr; + struct MaterialCylinder { + SurfaceId sid_; + CylPtr surface_; + std::string material_; + double thickness_; + MaterialCylinder(SurfaceId const& sid, CylPtr const& surface, std::string const& material, double thickness) : + sid_(sid), surface_(surface), material_(material), thickness_(thickness) {} + }; + using MaterialCylinderCollection = std::vector; // default constructor with nominal geometry DetectorSolenoid( CylPtr inner, CylPtr outer, DiskPtr front, DiskPtr back, CylPtr ipa, DiskPtr ipafront, DiskPtr ipaback, - FruPtr opa, AnnPtr tsda) : + FruPtr opa, AnnPtr tsda, MaterialCylinderCollection const& materialCylinders = MaterialCylinderCollection()) : inner_(inner) , outer_(outer), front_(front), back_(back), - ipa_(ipa), ipa_front_(ipafront), ipa_back_(ipaback), opa_(opa), tsda_(tsda) + ipa_(ipa), ipa_front_(ipafront), ipa_back_(ipaback), opa_(opa), tsda_(tsda), + materialCylinders_(materialCylinders) {} // accessors @@ -47,6 +59,7 @@ namespace mu2e { auto const& innerProtonAbsorberBackPtr() const { return ipa_back_; } auto const& outerProtonAbsorberPtr() const { return opa_; } auto const& upstreamAbsorberPtr() const { return tsda_; } + auto const& materialCylinders() const { return materialCylinders_; } private: CylPtr inner_; // inner cryostat cylinder boundary CylPtr outer_; // outer cryostat cylinder boundary @@ -57,6 +70,7 @@ namespace mu2e { DiskPtr ipa_back_; FruPtr opa_; // outer proton absorber AnnPtr tsda_; // TS downstream absorber + MaterialCylinderCollection materialCylinders_; // passive cylindrical material shells }; } } diff --git a/KinKalGeom/inc/KKMaterial.hh b/KinKalGeom/inc/KKMaterial.hh index f670b148f0..24a952c5ac 100644 --- a/KinKalGeom/inc/KKMaterial.hh +++ b/KinKalGeom/inc/KKMaterial.hh @@ -45,6 +45,12 @@ namespace mu2e { auto IPAMaterial() const { return matdbinfo_->findDetMaterial(ipamatname_); } auto STMaterial() const { return matdbinfo_->findDetMaterial(stmatname_); } auto CRVMaterial() const { return matdbinfo_->findDetMaterial(crvmatname_); } + auto material(std::string const& name) const { return matdbinfo_->findDetMaterial(name); } + // True when KinKal's ionization energy loss is the (restricted) Moyal mean, the only mode for + // which the KKShellXing unrestricted-Bethe-mean path correction is valid. Set from the + // IonizationEnergyLossMode fcl parameter so the correction stays consistent with the KinKal + // energy-loss model. See Mu2eKinKal/inc/KKShellXing.hh. + bool applyBetheCorrection() const { return betheCorrection_; } // FileFinder interface std::string matElmDictionaryFileName() const override; @@ -60,6 +66,7 @@ namespace mu2e { std::string wallmatname_, gasmatname_, wirematname_,ipamatname_, stmatname_, crvmatname_; std::unique_ptr matdbinfo_; // material database std::unique_ptr smat_; // straw material + bool betheCorrection_ = false; // KinKal eloss mode == Moyal mean (gates the KKShellXing Bethe correction) }; } #endif diff --git a/KinKalGeom/inc/KinKalGeom.hh b/KinKalGeom/inc/KinKalGeom.hh index 0c73a3b16f..3216e28652 100644 --- a/KinKalGeom/inc/KinKalGeom.hh +++ b/KinKalGeom/inc/KinKalGeom.hh @@ -11,18 +11,32 @@ #include "Offline/KinKalGeom/inc/CRV.hh" #include "Offline/DataProducts/inc/SurfaceId.hh" #include "KinKal/Geometry/Surface.hh" +#include "KinKal/Geometry/Plane.hh" +#include "KinKal/Geometry/Rectangle.hh" #include "Offline/Mu2eInterfaces/inc/Detector.hh" #include #include #include +#include namespace mu2e { class KinKalGeom : public Detector { public: using SurfacePtr = std::shared_ptr; + using RecPtr = std::shared_ptr; + using PlanePtr = std::shared_ptr; using SurfacePair =std::pair; using SurfacePairCollection = std::vector; using SurfacePairIter = std::multimap::const_iterator; using KKGMap = std::multimap; + struct PassiveMaterialPlane { + SurfaceId sid_; + PlanePtr surface_; // any planar surface: Rectangle (concrete/strongback/CRV) or Annulus (DS end walls) + std::string material_; + double thickness_; + PassiveMaterialPlane(SurfaceId const& sid, PlanePtr const& surface, std::string const& material, double thickness) : + sid_(sid), surface_(surface), material_(material), thickness_(thickness) {} + }; + using PassiveMaterialPlaneCollection = std::vector; // default constructor, now using GeometryService to create content KinKalGeom(){} // accessor to the raw map @@ -36,12 +50,14 @@ namespace mu2e { auto const& ST() const {return st_; } auto const& tracker() const {return tracker_; } auto const& CRV() const {return crv_; } + auto const& passiveMaterialPlanes() const { return passiveMaterialPlanes_; } private: // local copy of detector objects; these hold the actual (typed) surface objects std::unique_ptr tracker_; std::unique_ptr ds_; std::unique_ptr st_; std::unique_ptr crv_; + PassiveMaterialPlaneCollection passiveMaterialPlanes_; // the map used to find surfaces by Id KKGMap map_; // allow GeometryService to access internals, to avoid link loop diff --git a/KinKalGeom/src/KKMaterial.cc b/KinKalGeom/src/KKMaterial.cc index 0a88442830..ee8e0fd3d8 100644 --- a/KinKalGeom/src/KKMaterial.cc +++ b/KinKalGeom/src/KKMaterial.cc @@ -1,6 +1,7 @@ #include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/ConfigTools/inc/ConfigFileLookupPolicy.hh" #include "KinKal/MatEnv/DetMaterial.hh" +#include "messagefacility/MessageLogger/MessageLogger.h" namespace mu2e { using MatDBInfo = MatEnv::MatDBInfo; @@ -18,6 +19,12 @@ namespace mu2e { crvmatname_(matconfig.CRVMaterialName()) { MatEnv::DetMaterialConfig dmconf; dmconf.elossmode_ = (DetMaterial::energylossmode)matconfig.elossMode(); + // The KKShellXing Bethe (unrestricted ionization mean) correction is derived assuming KinKal + // returns the restricted Moyal mean; only enable it in that mode so the two stay consistent. + betheCorrection_ = (dmconf.elossmode_ == DetMaterial::moyalmean); + if(!betheCorrection_) + mf::LogInfo("KKMaterial") << "IonizationEnergyLossMode is not moyalmean (" + << matconfig.elossMode() << "): the KKShellXing Bethe-mean path correction is inactive."; dmconf.scatterfrac_solid_ = matconfig.solidScatter(); dmconf.scatterfrac_gas_ = matconfig.gasScatter(); dmconf.ebrehmsfrac_ = matconfig.eBrehms(); diff --git a/KinKalGeom/src/SConscript b/KinKalGeom/src/SConscript index 325b8e8f5c..eea1886bbb 100644 --- a/KinKalGeom/src/SConscript +++ b/KinKalGeom/src/SConscript @@ -23,6 +23,7 @@ mainlib = helper.make_mainlib([ 'MathCore', 'art_Framework_Services_Registry', 'art_Utilities', + 'MF_MessageLogger', 'canvas', 'cetlib', 'cetlib_except', diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index 90e691ed09..c4f32c80d4 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -170,7 +170,7 @@ namespace mu2e { // we have a good intersection. Use this to create a Shell material Xing auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); auto const& IPA = kkg_h->DS()->innerProtonAbsorberPtr(); - KKIPAXINGPTR ipaxingptr = std::make_shared(IPA,IPASID,*kkmat_h->IPAMaterial(),extrapIPA.intersection(),reftrajptr,ipathick_,extrapIPA.interTolerance()); + KKIPAXINGPTR ipaxingptr = std::make_shared(IPA,IPASID,*kkmat_h->IPAMaterial(),extrapIPA.intersection(),reftrajptr,ipathick_,extrapIPA.interTolerance(),kkmat_h->applyBetheCorrection()); if(extrapIPA.debug() > 2){ double dmom, paramomvar, perpmomvar; ipaxingptr->materialEffects(dmom,paramomvar,perpmomvar); @@ -211,7 +211,7 @@ namespace mu2e { if(extrapST.intersection().good()){ // we have a good intersection. Use this to create a Shell material Xing auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); - KKSTXINGPTR stxingptr = std::make_shared(extrapST.foil(),extrapST.foilId(),*kkmat_h->STMaterial(),extrapST.intersection(),reftrajptr,stthick_,extrapST.interTolerance()); + KKSTXINGPTR stxingptr = std::make_shared(extrapST.foil(),extrapST.foilId(),*kkmat_h->STMaterial(),extrapST.intersection(),reftrajptr,stthick_,extrapST.interTolerance(),kkmat_h->applyBetheCorrection()); if(extrapST.debug() > 2){ double dmom, paramomvar, perpmomvar; stxingptr->materialEffects(dmom,paramomvar,perpmomvar); @@ -306,7 +306,7 @@ namespace mu2e { auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); auto crvxingptr = std::make_shared(extrapCRV.sector((size_t)inter.isect_).sector_, SurfaceId(kkg_h->CRV()->sectorName(inter.isect_)),*kkmat_h->CRVMaterial(),inter.inter_,reftrajptr, - 2*inter.whw_, extrapCRV.interTolerance()); + 2*inter.whw_, extrapCRV.interTolerance(),kkmat_h->applyBetheCorrection()); ktrk.addCRVXing(crvxingptr,tdir); if(debug_ > 1) std::cout << "Good CRV " << inter.inter_ << ftraj.range() << std::endl; } diff --git a/Mu2eKinKal/inc/KKShellXing.hh b/Mu2eKinKal/inc/KKShellXing.hh index c53f3c2294..83420fb113 100644 --- a/Mu2eKinKal/inc/KKShellXing.hh +++ b/Mu2eKinKal/inc/KKShellXing.hh @@ -9,8 +9,40 @@ #include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" #include "KinKal/MatEnv/DetMaterial.hh" #include "Offline/DataProducts/inc/SurfaceId.hh" +#include +#include namespace mu2e { + // Ratio of the unrestricted Bethe ionization mean to KinKal's energyLoss (the restricted/moyal loss), >=1. + // We correct the eloss for thick passive materials (ie. concrete blocks) by scaling the crossing path length. + // Side-effect: the multiple-scattering variance scales too (~3-8% larger pointing sigma at high p), + // which does not touch the |p| residual. f>=1 always (we never reduce the loss). + // Only applied if KinKal is configured with IonizationEnergyLossMode = moyalmean (KKMaterial fcl). + // + // This hand-rolls the unrestricted Bethe mean from DetMaterial's eloss_xi/eexc/densityCorrection. + // It intentionally omits DetMaterial::shellCorrection (which ionizationEnergyLossMPV includes): that + // term is negligible for the ~100 MeV/c muons and electrons extrapolated here, and this whole helper + // is a stopgap. + // TODO(remove once KinKal provides the unrestricted ionization mean directly as an energy-loss option, + // then this ratio and the betheCorr plumbing can be deleted). + inline double betheCorrectionFactor(MatEnv::DetMaterial const& mat, double mom, double pathlen, double mass) { + if(mom <= 0.0 || pathlen == 0.0) return 1.0; + static constexpr double me = 0.510998950; // electron mass [MeV] + double const beta2 = mom*mom/(mom*mom + mass*mass); + double const bg2 = (mom*mom)/(mass*mass); // (beta*gamma)^2 + double const gamma = std::sqrt(mom*mom + mass*mass)/mass; + double const xi = mat.eloss_xi(std::sqrt(beta2), pathlen); // (K/2)(Z/A)*rho*|path|/beta^2 [MeV] + double const I = mat.eexc(); // mean excitation energy [MeV] + double const Tmax = 2.0*me*bg2/(1.0 + 2.0*gamma*me/mass + (me/mass)*(me/mass)); + double const delta = mat.densityCorrection(bg2); // same density-effect KinKal uses + // unrestricted Bethe MEAN loss, as a NEGATIVE energy change (DetMaterial::energyLoss is negative for loss) + double const meanChange = -xi*( std::log(2.0*me*bg2/I) + std::log(Tmax/I) - 2.0*beta2 - delta ); + double const el = mat.energyLoss(mom, pathlen, mass); // moyalmean, negative + if(el >= 0.0) return 1.0; + double const f = meanChange/el; // |mean|/|moyal|, both negative -> positive + return f > 1.0 ? f : 1.0; // only ever increase the loss + } + template class KKShellXing : public KinKal::ElementXing { public: using PTRAJ = KinKal::ParticleTrajectory; @@ -19,8 +51,8 @@ namespace mu2e { using PCA = KinKal::PiecewiseClosestApproach; using CA = KinKal::ClosestApproach; using SURFPTR = std::shared_ptr; - // construct from a surface, material, intersection, and transverse thickness - KKShellXing(SURFPTR surface, SurfaceId const& sid, MatEnv::DetMaterial const& mat, KinKal::Intersection inter, KTRAJPTR reftraj, double thickness, double tol); + // construct from a surface, material, intersection, and transverse thickness. + KKShellXing(SURFPTR surface, SurfaceId const& sid, MatEnv::DetMaterial const& mat, KinKal::Intersection inter, KTRAJPTR reftraj, double thickness, double tol, bool betheCorr); virtual ~KKShellXing() {} // clone op for reinstantiation KKShellXing(KKShellXing const& rhs) = default; @@ -56,19 +88,24 @@ namespace mu2e { std::vector mxings_; // material xing double thick_; // shell thickness double tol_; // tolerance for intersection + bool betheCorr_; // apply the unrestricted-Bethe-mean eloss path correction (Moyal-mean mode only) double varscale_; // variance scale, for annealing KinKal::Parameters fparams_; // 1st-order parameter change for forwards time }; template KKShellXing::KKShellXing(SURFPTR surface, SurfaceId const& sid, MatEnv::DetMaterial const& mat, KinKal::Intersection inter, - std::shared_ptr reftrajptr, double thickness, double tol) : + std::shared_ptr reftrajptr, double thickness, double tol, bool betheCorr) : surf_(surface), sid_(sid), mat_(mat), inter_(inter), reftrajptr_(reftrajptr), thick_(thickness),tol_(tol), - varscale_(1.0) + betheCorr_(betheCorr), varscale_(1.0) { if(inter_.good()){ // compute the path length - double pathlen = thick_/(inter_.norm_.Dot(inter_.pdir_)); - mxings_.emplace_back(mat_,pathlen); + double dotprod = std::max(1e-6,std::fabs(inter_.norm_.Dot(inter_.pdir_))); + double pathlen = thick_/dotprod; + // Scale the crossed path length so the (restricted/moyal) E loss becomes the unrestricted Bethe mean. + // Only when the KinKal eloss model is the Moyal mean (betheCorr_); otherwise leave the loss unscaled. + double const f = betheCorr_ ? betheCorrectionFactor(mat_, reftrajptr_->momentum(), pathlen, reftrajptr_->mass()) : 1.0; + mxings_.emplace_back(mat_, pathlen*f); } } @@ -90,8 +127,11 @@ namespace mu2e { // check if we are on the surface; if so, create the xing if(inter_.good()){ // compute the path length - double pathlen = thick_/(inter_.norm_.Dot(inter_.pdir_)); - mxings_.emplace_back(mat_,pathlen); + double dotprod = std::max(1e-6,std::fabs(inter_.norm_.Dot(inter_.pdir_))); + double pathlen = thick_/dotprod; + // same Bethe-mean path scale as the ctor (keeps the fit-side params consistent with the extrapolation) + double const f = betheCorr_ ? betheCorrectionFactor(mat_, reftrajptr_->momentum(), pathlen, reftrajptr_->mass()) : 1.0; + mxings_.emplace_back(mat_, pathlen*f); fparams_ = this->parameterChange(varscale_); } } diff --git a/TrackerConditions/data/MaterialsList.data b/TrackerConditions/data/MaterialsList.data index 2d64c9a003..52f9b9fd01 100644 --- a/TrackerConditions/data/MaterialsList.data +++ b/TrackerConditions/data/MaterialsList.data @@ -30,15 +30,34 @@ cathode 5.2 0. 0. +2 86.4e-2 Copper 0 13.6e-2 straw-wall 1 -10 straw_15um 1.927e-2 0.0 0.0 +2 91.3e-2 straw-wall 1 8.7e-2 straw-gas 1 -10 -20 -30 20.0 1.0 solid straw_25um 3.068e-2 0.0 0.0 +2 94.6e-2 straw-wall 1 5.4e-2 straw-gas 1 -10 -20 -30 20.0 1.0 solid # DS materials -NbTi 6.5 0.0 0.0 +2 0.65 Niobium 0.35 Titanium 0 -10 -20 -30 20.0 1.0 solid -DS1Coil 3.35 0.0 0.0 +2 0.67 Aluminum 0.33 NbTi 0 -10 -20 -30 20.0 1.0 solid -DS2Coil 3.031 0.0 0.0 +2 0.813 Aluminum 0.187 NbTi 0 -10 -20 -30 20.0 1.0 solid -DSSStainless 8.02 0.0 0.0 +5 0.02 Manganese 0.01 Silicon 0.19 Chromium 0.1 Nickel 0.68 Iron 0 -10 -20 -30 20.0 1.0 solid +NbTi 6.5 0.0 0.0 +2 0.65 Niobium 0 0.35 Titanium 0 -10 -20 -30 20.0 1.0 solid +# Pure aluminium (2.700 g/cm3) passive materials, named per subsystem so the surface/material is +# self-documenting in the KinKalGeom passive-plane dumps. Same composition as the stopping-target +# 'Target' material; kept distinct for clarity. +CRVStrongbackAl 2.700000 0.0 0.0 +1 100.0e-2 Aluminum 0 24.01000 106.40000 -30 20.0 1.0 solid +DSShieldAl 2.700000 0.0 0.0 +1 100.0e-2 Aluminum 0 24.01000 106.40000 -30 20.0 1.0 solid +# DS1Coil / DS2Coil: individual coil materials, kept here for reference only - the extrapolation uses +# the collapsed DSCoilMix blend below (see its derivation), not these directly. +#DS1Coil 3.35 0.0 0.0 +2 0.67 Aluminum 0 0.33 NbTi 1 -10 -20 -30 20.0 1.0 solid +#DS2Coil 3.031 0.0 0.0 +2 0.813 Aluminum 0 0.187 NbTi 1 -10 -20 -30 20.0 1.0 solid +# DSCoilMix: the 11 DS coils collapsed into ONE effective material for the KinKalGeom long-cylinder +# coil approximation (individual coils are too short in z to be reliably intersected by the extrapolation). +# z-coverage-weighted blend of the 8 DS1CoilMix coils (41 mm, 3.35 g/cm3, 67/33 Al/NbTi, z-cover 3128 mm) +# and the 3 DS2CoilMix coils (20.5 mm, 3.031 g/cm3, 81.3/18.7 Al/NbTi, z-cover 5332.5 mm): z-mass composition +# 73.23% Al + 26.77% NbTi, areal mass 8.994 g/cm2; density = areal mass / (z-weighted thickness 28.08 mm) = +# 3.2032, scaled to the geometric thickness as CRVModule is. Reproduces the true z-avg t/X0 = 0.473. +DSCoilMix 3.203230 0.0 0.0 +2 0.73227 Aluminum 0 0.26773 NbTi 1 -10 -20 -30 20.0 1.0 solid +StainlessSteel 8.02 0.0 0.0 +5 0.02 Manganese 0 0.01 Silicon 0 0.19 Chromium 0 0.1 Nickel 0 0.68 Iron 0 -10 -20 -30 20.0 1.0 solid +# Concrete definitions translated from Mu2eG4/src/ConstructMaterials.cc. +CONCRETE_MARS 2.35 0.0 0.0 +9 0.006 Hydrogen 0 0.030 Carbon 0 0.500 Oxygen 0 0.010 Sodium 0 0.030 Aluminum 0 0.200 Silicon 0 0.010 Potasium 0 0.200 Calcium 0 0.014 Iron 0 -10 -20 -30 20.0 1.0 solid +CONCRETE_MASONRY 1.16 0.0 0.0 +9 0.006 Hydrogen 0 0.030 Carbon 0 0.500 Oxygen 0 0.010 Sodium 0 0.030 Aluminum 0 0.200 Silicon 0 0.010 Potasium 0 0.200 Calcium 0 0.014 Iron 0 -10 -20 -30 20.0 1.0 solid # CRV active module stack from crv_parameters.xlsx (DocDB 862 https://mu2e-docdb.fnal.gov/cgi-bin/sso/ShowDocument?docid=862): # 4 scintillator layers (4*19.80 mm) + 3 aluminum absorbers (3*9.525 mm) -# + aluminum strongback (12.7 mm) + thin aluminum cover (3.175 mm), +# + thin aluminum cover (3.175 mm). The aluminum strongback (12.7 mm) is modeled +# separately as the CRV_StrongBack passive plane, so it is excluded here to avoid +# double-counting its mass and scattering. # represented as one effective material over the KinKal CRV thickness. -# Multiply density by (full module width)/(KinKal CRV thickness) = 123.65 / 107.775 +# Multiply density by (module width excl. strongback)/(KinKal CRV thickness) = 110.95 / 107.775 # CRVModule density is computed from this areal density and rounded to 6 decimal places. Polystyrene 1.060000 0. 0. -2 8 Carbon 0 8 Hydrogen 0 -10 -20 -30 20.0 1.0 solid -CRVModule 1.892526 0. 0. +2 41.2e-2 Polystyrene 1 58.8e-2 Aluminum 0 -10 -20 -30 20.0 1.0 solid +CRVModule 1.574363 0. 0. +2 49.4775e-2 Polystyrene 1 50.5225e-2 Aluminum 0 -10 -20 -30 20.0 1.0 solid