From f2f88991d2c72dfe3e778d96141c6e2db7d78b68 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Mon, 11 May 2026 16:01:26 -0700 Subject: [PATCH 01/12] Add materials for DS --- GeometryService/src/KinKalGeomMaker.cc | 10 ++++++++++ TrackerConditions/data/MaterialsList.data | 7 ++++++- 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/GeometryService/src/KinKalGeomMaker.cc b/GeometryService/src/KinKalGeomMaker.cc index 5c955ef72e..e05f1c128a 100644 --- a/GeometryService/src/KinKalGeomMaker.cc +++ b/GeometryService/src/KinKalGeomMaker.cc @@ -79,6 +79,16 @@ namespace mu2e { } void KinKalGeomMaker::makeDS() { + GeomHandle det; + GeomHandle ds; + std::cout << "DS Cryo or " << ds->rOut1() << "," << ds->rOut2() << " ir " << ds->rIn1()<<","<< ds->rIn2() << " halfl " << ds->halfLength() + << " zpos " << ds->position().z() << " material " << ds->material() << std::endl; + std::cout << "DS shield or " << ds->shield_rOut1() << "," << ds->shield_rOut2() << " ir " << ds->shield_rIn1()<<","<< ds->shield_rIn2() << " halfl " << ds->shield_halfLength() << " material " << ds->shield_material() << std::endl; + std::cout << "DS ncoils " << ds->nCoils() << std::endl; + for(size_t icoil = 0; icoil < static_cast(ds->nCoils()); icoil++){ + std::cout << "DS coil ir " << ds->coil_rIn() << " or " << ds->coil_rOut()[icoil] << " length " << ds->coil_zLength()[icoil] << " zpos " << ds->coil_zPosition()[icoil] + << " material " << ds->coil_materials()[icoil] << std::endl; + } // currently use hard-coded geometry auto inner= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-1482),950,5450); auto outer= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-1482),1328,5450); // bounding surfaces diff --git a/TrackerConditions/data/MaterialsList.data b/TrackerConditions/data/MaterialsList.data index 9d2b64e7e6..01b9b5451b 100644 --- a/TrackerConditions/data/MaterialsList.data +++ b/TrackerConditions/data/MaterialsList.data @@ -1,5 +1,5 @@ # -# This is an exhaustive list of materials for entry into the Conditions DB +# This is an exhaustive list of materials used in Mu2e KinKal # # Do not remove the following line it's a tag ! #Materials_list @@ -29,4 +29,9 @@ cathode 5.2 0. 0. +2 86.4e-2 Copper 0 13.6e-2 straw-wall 1 -10 # effective straw wall + gas material, assuming uniform mass distribution 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 From bfc1fd998fd24e885f6b0527c38016350134af7b Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Sat, 16 May 2026 10:32:20 -0700 Subject: [PATCH 02/12] Prep for moving KKMaterial to KinKalGeom --- GeometryService/src/KinKalGeomMaker.cc | 26 +++++++++++++++++++++++--- KinKalGeom/inc/DetectorSolenoid.hh | 4 ++-- KinKalGeom/inc/KinKalGeom.hh | 2 +- KinKalGeom/inc/ShellMaterial.hh | 21 +++++++++++++++++++++ Mu2eKinKal/inc/KKStrawMaterial.hh | 10 ++++------ Mu2eKinKal/inc/KKStrawXing.hh | 5 +++-- Mu2eKinKal/src/KKStrawMaterial.cc | 16 +++++++--------- RecoDataProducts/inc/TrkStraw.hh | 3 ++- 8 files changed, 63 insertions(+), 24 deletions(-) create mode 100644 KinKalGeom/inc/ShellMaterial.hh diff --git a/GeometryService/src/KinKalGeomMaker.cc b/GeometryService/src/KinKalGeomMaker.cc index e05f1c128a..9d33dbccdb 100644 --- a/GeometryService/src/KinKalGeomMaker.cc +++ b/GeometryService/src/KinKalGeomMaker.cc @@ -89,11 +89,31 @@ namespace mu2e { std::cout << "DS coil ir " << ds->coil_rIn() << " or " << ds->coil_rOut()[icoil] << " length " << ds->coil_zLength()[icoil] << " zpos " << ds->coil_zPosition()[icoil] << " material " << ds->coil_materials()[icoil] << std::endl; } - // currently use hard-coded geometry - auto inner= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-1482),950,5450); - auto outer= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-1482),1328,5450); // bounding surfaces + //DS Cryo or 1303,1328 ir 950,969.05 halfl 5450 zpos 8689 material StainlessSteel + //DS shield or 1237.3,1250 ir 1010,1022.7 halfl 5287.7 material G4_Al + //DS ncoils 11 + //DS coil ir 1050 or 1091 length 419.75 zpos 3539 material DS1CoilMix + //DS coil ir 1050 or 1091 length 419.75 zpos 3964 material DS1CoilMix + //DS coil ir 1050 or 1091 length 419.75 zpos 4389 material DS1CoilMix + //DS coil ir 1050 or 1091 length 419.75 zpos 5042 material DS1CoilMix + //DS coil ir 1050 or 1091 length 362.25 zpos 5699 material DS1CoilMix + //DS coil ir 1050 or 1091 length 362.25 zpos 6369 material DS1CoilMix + //DS coil ir 1050 or 1091 length 362.25 zpos 7176 material DS1CoilMix + //DS coil ir 1050 or 1070.5 length 1777.5 zpos 7949 material DS2CoilMix + //DS coil ir 1050 or 1070.5 length 1777.5 zpos 9761 material DS2CoilMix + //DS coil ir 1050 or 1070.5 length 1777.5 zpos 11544 material DS2CoilMix + //DS coil ir 1050 or 1091 length 362.25 zpos 13326 material DS1CoilMix + // + // + // + // + auto inner= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-1482),ds->rIn1(),ds->halfLength()); + 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()); + + + // 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); auto ipafront= std::make_shared(ipa->frontDisk()); auto ipaback= std::make_shared(ipa->backDisk()); diff --git a/KinKalGeom/inc/DetectorSolenoid.hh b/KinKalGeom/inc/DetectorSolenoid.hh index 7721710689..5fd98a55db 100644 --- a/KinKalGeom/inc/DetectorSolenoid.hh +++ b/KinKalGeom/inc/DetectorSolenoid.hh @@ -48,8 +48,8 @@ namespace mu2e { auto const& outerProtonAbsorberPtr() const { return opa_; } auto const& upstreamAbsorberPtr() const { return tsda_; } private: - CylPtr inner_; // inner cryostat cylinder - CylPtr outer_; // outer cryostat cylinder + CylPtr inner_; // inner cryostat cylinder boundary + CylPtr outer_; // outer cryostat cylinder boundary DiskPtr front_; // front (upstream) and back (downstream) of DS DiskPtr back_; CylPtr ipa_; // inner proton absorber diff --git a/KinKalGeom/inc/KinKalGeom.hh b/KinKalGeom/inc/KinKalGeom.hh index 05770e6730..3d9a7e7dc2 100644 --- a/KinKalGeom/inc/KinKalGeom.hh +++ b/KinKalGeom/inc/KinKalGeom.hh @@ -25,7 +25,7 @@ namespace mu2e { using SurfacePairCollection = std::vector; using SurfacePairIter = std::multimap::const_iterator; using KKGMap = std::multimap; - // default constructor, now using GeometryService + // default constructor, now using GeometryService to create content KinKalGeom(); // accessor to the raw map auto const& map() const { return map_; } diff --git a/KinKalGeom/inc/ShellMaterial.hh b/KinKalGeom/inc/ShellMaterial.hh new file mode 100644 index 0000000000..91a8338949 --- /dev/null +++ b/KinKalGeom/inc/ShellMaterial.hh @@ -0,0 +1,21 @@ +// +// Class modeling a thin material as a surface, material, and thickness. Can be used to create a ShellXing in KKTrack +// Original author: David Brown (LBNL) 5/26 +// +#ifndef KinKalGeom_ShellMaterial_hh +#define KinKalGeom_ShellMaterial_hh +#include "Offline/DataProducts/inc/SurfaceId.hh" +#include "KinKal/Geometry/Surface.hh" +#include "KinKal/MatEnv/DetMaterial.hh" + +namespace mu2e { + class ShellMaterial { + using SurfacePtr = std::shared_ptr; + public: + private: + + SurfacePtr surf_; // surface for this materia + + + }; +} diff --git a/Mu2eKinKal/inc/KKStrawMaterial.hh b/Mu2eKinKal/inc/KKStrawMaterial.hh index ab76fe02b1..d6cdc18381 100644 --- a/Mu2eKinKal/inc/KKStrawMaterial.hh +++ b/Mu2eKinKal/inc/KKStrawMaterial.hh @@ -6,7 +6,6 @@ // #include "KinKal/Trajectory/ClosestApproachData.hh" #include "Offline/TrackerGeom/inc/StrawProperties.hh" -#include "Offline/Mu2eKinKal/inc/StrawXingUpdater.hh" namespace MatEnv { class MatDBInfo; @@ -30,15 +29,14 @@ namespace mu2e { const std::shared_ptr gasmat_, const std::shared_ptr wiremat_); // construct using materials by name - KKStrawMaterial(MatDBInfo const& matdbinfo,StrawProperties const& sprops, - const std::string& wallmat="straw-wall", const std::string& gasmat="straw-gas", const std::string& wiremat="straw-wire"); + KKStrawMaterial(MatDBInfo const& matdbinfo,StrawProperties const& sprops); // pathlength through straw components, given closest approach. Return the method used to compute the paths - PathCalc pathLengths(ClosestApproachData const& cadata,StrawXingUpdater const& caconfig, double& wallpath, double& gaspath, double& wirepath) const; + PathCalc pathLengths(ClosestApproachData const& cadata,double nsig, double& wallpath, double& gaspath, double& wirepath,int diag=0) const; PathCalc averagePathLengths(double& wallpath, double& gaspath, double& wirepath) const; // transit length given closest approach double transitLength(ClosestApproachData const& cadata) const; - // find the material crossings given doca and error on doca. Should allow for straw and wire to have different axes TODO - PathCalc findXings(ClosestApproachData const& cadata,StrawXingUpdater const& caconfig, std::vector& mxings) const; + // find the material crossings given doca and error on doca. The CA object should be WRT the straw axis, not the wire axis + PathCalc findXings(ClosestApproachData const& cadata, double nsig, std::vector& mxings, int diag=0) const; DetMaterial const& wallMaterial() const { return *wallmat_; } DetMaterial const& gasMaterial() const { return *gasmat_; } DetMaterial const& wireMaterial() const { return *wiremat_; } diff --git a/Mu2eKinKal/inc/KKStrawXing.hh b/Mu2eKinKal/inc/KKStrawXing.hh index 4c09250ea3..e32f7e0239 100644 --- a/Mu2eKinKal/inc/KKStrawXing.hh +++ b/Mu2eKinKal/inc/KKStrawXing.hh @@ -7,10 +7,11 @@ #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/StrawXingUpdater.hh" #include "Offline/Mu2eKinKal/inc/KKStrawMaterial.hh" +#include "Offline/TrackerGeom/inc/Straw.hh" +#include "Offline/DataProducts/inc/StrawId.hh" #include "KinKal/Trajectory/ParticleTrajectory.hh" #include "KinKal/Trajectory/SensorLine.hh" #include "KinKal/Trajectory/PiecewiseClosestApproach.hh" -#include "Offline/DataProducts/inc/StrawId.hh" #include "cetlib_except/exception.h" namespace mu2e { using KinKal::SVEC3; @@ -164,7 +165,7 @@ namespace mu2e { varscale_ = 1.0; } // update the material xings from gas, straw wall, and wire - smat_.findXings(ca_.tpData(),sxconfig_,mxings_); + smat_.findXings(ca_.tpData(),sxconfig_.nsig_,mxings_,sxconfig_.diag_); // update the effect these have on the parameters fparams_ = this->parameterChange(varscale_); } diff --git a/Mu2eKinKal/src/KKStrawMaterial.cc b/Mu2eKinKal/src/KKStrawMaterial.cc index a79a54cef3..404ff4220d 100644 --- a/Mu2eKinKal/src/KKStrawMaterial.cc +++ b/Mu2eKinKal/src/KKStrawMaterial.cc @@ -35,11 +35,11 @@ namespace mu2e { return KKStrawMaterial::average; } - KKStrawMaterial::PathCalc KKStrawMaterial::pathLengths(ClosestApproachData const& cadata,StrawXingUpdater const& caconfig, - double& wallpath, double& gaspath, double& wirepath) const { + KKStrawMaterial::PathCalc KKStrawMaterial::pathLengths(ClosestApproachData const& cadata,double nsig, + double& wallpath, double& gaspath, double& wirepath,int diag) const { wallpath = gaspath = wirepath = 0.0; PathCalc retval = KKStrawMaterial::unknown; - double docarange = caconfig.nsig_*sqrt(std::max(0.0,cadata.docaVar())); + double docarange = nsig*sqrt(std::max(0.0,cadata.docaVar())); // if the doca range covers the straw, use the average double adoca = fabs(cadata.doca()); double mindoca = std::max(0.0,std::min(adoca,irad_)-docarange); @@ -58,10 +58,9 @@ namespace mu2e { if(awall>0.0)wallpath = awall/(omaxdoca-mindoca); retval = range; } - if(caconfig.diag_>0){ + if(diag>0){ std::cout << "KKStrawMaterial: DOCA " << fabs(cadata.doca()) << " DOCA Var " << cadata.docaVar() << " range [" << mindoca << "," << imaxdoca << "] wall path (mm)" << wallpath << " gas path " << gaspath << " retval " << retval << std::endl; } - // Model the wire as a diffuse gas, density constrained by DOCA TODO // 3D pathlength includes projection along straw double afac = angleFactor(cadata.dirDot()); @@ -79,15 +78,14 @@ namespace mu2e { return tlen; } - KKStrawMaterial::PathCalc KKStrawMaterial::findXings(ClosestApproachData const& cadata,StrawXingUpdater const& caconfig, - std::vector& mxings) const { + KKStrawMaterial::PathCalc KKStrawMaterial::findXings(ClosestApproachData const& cadata, double nsig, std::vector 0.0) mxings.push_back(MaterialXing(*wallmat_,wallpath)); if(gaspath > 0.0) mxings.push_back(MaterialXing(*gasmat_,gaspath)); if(wirepath > 0.0) mxings.push_back(MaterialXing(*wiremat_,wirepath)); - if(caconfig.diag_ > 1){ + if(diag > 1){ std::cout << "Ce wall dE/dx " << wallmat_->energyLoss(105,1.0,0.511) << std::endl; std::cout << "Ce gas dE/dx " << gasmat_->energyLoss(105,1.0,0.511) << std::endl; } diff --git a/RecoDataProducts/inc/TrkStraw.hh b/RecoDataProducts/inc/TrkStraw.hh index 3af028030b..59d4a6db99 100644 --- a/RecoDataProducts/inc/TrkStraw.hh +++ b/RecoDataProducts/inc/TrkStraw.hh @@ -9,6 +9,7 @@ #include "Offline/DataProducts/inc/StrawId.hh" #include "Offline/RecoDataProducts/inc/StrawFlag.hh" #include "Offline/Mu2eKinKal/inc/KKStrawMaterial.hh" +#include "Offline/Mu2eKinKal/inc/StrawXingUpdater.hh" #include "KinKal/Trajectory/ClosestApproachData.hh" namespace mu2e { @@ -25,7 +26,7 @@ namespace mu2e { _dmom(dmom) { double wallpath, gaspath, wirepath; - _pcalc = smat.pathLengths(pocadata,caconfig,wallpath,gaspath,wirepath); + _pcalc = smat.pathLengths(pocadata,caconfig.nsig_,wallpath,gaspath,wirepath); _wallpath = wallpath; _gaspath = gaspath; _wirepath = wirepath; From d792d4d8c8b5e485c079559ae971b7cf96b4d5e7 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Sat, 16 May 2026 10:52:41 -0700 Subject: [PATCH 03/12] Fix bugs --- Mu2eKinKal/inc/KKStrawMaterial.hh | 4 +++- Mu2eKinKal/src/KKStrawMaterial.cc | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/Mu2eKinKal/inc/KKStrawMaterial.hh b/Mu2eKinKal/inc/KKStrawMaterial.hh index d6cdc18381..d29f153a12 100644 --- a/Mu2eKinKal/inc/KKStrawMaterial.hh +++ b/Mu2eKinKal/inc/KKStrawMaterial.hh @@ -6,6 +6,7 @@ // #include "KinKal/Trajectory/ClosestApproachData.hh" #include "Offline/TrackerGeom/inc/StrawProperties.hh" +#include namespace MatEnv { class MatDBInfo; @@ -29,7 +30,8 @@ namespace mu2e { const std::shared_ptr gasmat_, const std::shared_ptr wiremat_); // construct using materials by name - KKStrawMaterial(MatDBInfo const& matdbinfo,StrawProperties const& sprops); + KKStrawMaterial(MatDBInfo const& matdbinfo,StrawProperties const& sprops, + const std::string& wallmat, const std::string& gasmat, const std::string& wiremat); // pathlength through straw components, given closest approach. Return the method used to compute the paths PathCalc pathLengths(ClosestApproachData const& cadata,double nsig, double& wallpath, double& gaspath, double& wirepath,int diag=0) const; PathCalc averagePathLengths(double& wallpath, double& gaspath, double& wirepath) const; diff --git a/Mu2eKinKal/src/KKStrawMaterial.cc b/Mu2eKinKal/src/KKStrawMaterial.cc index 404ff4220d..003d878fd3 100644 --- a/Mu2eKinKal/src/KKStrawMaterial.cc +++ b/Mu2eKinKal/src/KKStrawMaterial.cc @@ -78,7 +78,7 @@ namespace mu2e { return tlen; } - KKStrawMaterial::PathCalc KKStrawMaterial::findXings(ClosestApproachData const& cadata, double nsig, std::vector& mxings, int diag) const { mxings.clear(); double wallpath, gaspath, wirepath; auto retval = pathLengths(cadata,nsig,wallpath, gaspath, wirepath,diag); From 3c30b9ccc5730e7fde3189a99ffc08a72edeb65c Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Sat, 16 May 2026 11:31:38 -0700 Subject: [PATCH 04/12] Consolidate interfaces --- Mu2eKinKal/inc/KKFileFinder.hh | 48 ---------------------------------- Mu2eKinKal/inc/KKMaterial.hh | 22 +++++++++++++--- Mu2eKinKal/src/KKMaterial.cc | 13 ++++++--- 3 files changed, 28 insertions(+), 55 deletions(-) delete mode 100644 Mu2eKinKal/inc/KKFileFinder.hh diff --git a/Mu2eKinKal/inc/KKFileFinder.hh b/Mu2eKinKal/inc/KKFileFinder.hh deleted file mode 100644 index 67ab289602..0000000000 --- a/Mu2eKinKal/inc/KKFileFinder.hh +++ /dev/null @@ -1,48 +0,0 @@ -#ifndef Mu2eKinKal_KKFileFinder_hh -#define Mu2eKinKal_KKFileFinder_hh - -#include "KinKal/MatEnv/FileFinderInterface.hh" -#include "Offline/ConfigTools/inc/ConfigFileLookupPolicy.hh" -#include - -namespace mu2e { - - class KKFileFinder : public MatEnv::FileFinderInterface { - - public: - KKFileFinder(std::string elementsBaseName, - std::string isotopesBaseName, - std::string materialsBaseName): - elementsBaseName_(elementsBaseName), - isotopesBaseName_(isotopesBaseName), - materialsBaseName_(materialsBaseName){} - - std::string matElmDictionaryFileName() const override { - return findFile(elementsBaseName_); - } - - std::string matIsoDictionaryFileName() const override { - return findFile(isotopesBaseName_); - } - - std::string matMtrDictionaryFileName() const override { - return findFile(materialsBaseName_); - } - - std::string findFile( std::string const& path ) const override { - return policy_( path ); - } - - private: - - mutable ConfigFileLookupPolicy policy_; - - std::string elementsBaseName_; - std::string isotopesBaseName_; - std::string materialsBaseName_; - - }; - -} - -#endif /* btrkHelper_KKFileFinder_hh */ diff --git a/Mu2eKinKal/inc/KKMaterial.hh b/Mu2eKinKal/inc/KKMaterial.hh index 127d59642d..fecefe8e04 100644 --- a/Mu2eKinKal/inc/KKMaterial.hh +++ b/Mu2eKinKal/inc/KKMaterial.hh @@ -8,14 +8,17 @@ #include "fhiclcpp/types/Tuple.h" // KinKal #include "KinKal/MatEnv/MatDBInfo.hh" +#include "KinKal/MatEnv/FileFinderInterface.hh" +// Mu2eKinKal: moveme! #include "Offline/Mu2eKinKal/inc/KKStrawMaterial.hh" -#include "Offline/Mu2eKinKal/inc/KKFileFinder.hh" +// mu2e +#include "Offline/ConfigTools/inc/ConfigFileLookupPolicy.hh" #include #include namespace mu2e { - class KKMaterial { + class KKMaterial : public MatEnv::FileFinderInterface { public: using Name = fhicl::Name; using Comment = fhicl::Comment; @@ -40,11 +43,22 @@ namespace mu2e { KKStrawMaterial const& strawMaterial() const; auto IPAMaterial() const { return matdbinfo_->findDetMaterial(ipamatname_); } auto STMaterial() const { return matdbinfo_->findDetMaterial(stmatname_); } + + // FileFinder interface + std::string matElmDictionaryFileName() const override; + std::string matIsoDictionaryFileName() const override; + std::string matMtrDictionaryFileName() const override; + std::string findFile( std::string const& path ) const override; private: - KKFileFinder filefinder_; // used to find material info + // material description files base names (not path) + std::string elementsBaseName_; + std::string isotopesBaseName_; + std::string materialsBaseName_; + mutable ConfigFileLookupPolicy policy_; + // specific material names std::string wallmatname_, gasmatname_, wirematname_,ipamatname_, stmatname_; mutable std::unique_ptr matdbinfo_; // material database - mutable std::unique_ptr smat_; // straw material + mutable std::unique_ptr smat_; // straw material; move to KinKalGeom }; } #endif diff --git a/Mu2eKinKal/src/KKMaterial.cc b/Mu2eKinKal/src/KKMaterial.cc index 126193fc9b..72ff52ac6c 100644 --- a/Mu2eKinKal/src/KKMaterial.cc +++ b/Mu2eKinKal/src/KKMaterial.cc @@ -2,14 +2,15 @@ #include "Offline/GeometryService/inc/GeomHandle.hh" #include "Offline/TrackerGeom/inc/Tracker.hh" #include "KinKal/MatEnv/DetMaterial.hh" -#include "Offline/GeometryService/inc/GeometryService.hh" namespace mu2e { using MatDBInfo = MatEnv::MatDBInfo; using MatEnv::DetMaterial; KKMaterial::KKMaterial(KKMaterial::Config const& matconfig) : - filefinder_(matconfig.elements(),matconfig.isotopes(),matconfig.materials()), + elementsBaseName_(matconfig.elements()), + isotopesBaseName_(matconfig.isotopes()), + materialsBaseName_(matconfig.materials()), wallmatname_(matconfig.strawWallMaterialName()), gasmatname_(matconfig.strawGasMaterialName()), wirematname_(matconfig.strawWireMaterialName()), @@ -20,7 +21,7 @@ namespace mu2e { dmconf.scatterfrac_solid_ = matconfig.solidScatter(); dmconf.scatterfrac_gas_ = matconfig.gasScatter(); dmconf.ebrehmsfrac_ = matconfig.eBrehms(); - matdbinfo_ = std::make_unique(filefinder_,dmconf); + matdbinfo_ = std::make_unique(*this,dmconf); } KKStrawMaterial const& KKMaterial::strawMaterial() const { @@ -36,4 +37,10 @@ namespace mu2e { } return *smat_; } + + std::string KKMaterial::findFile( std::string const& basename ) const { return policy_( basename ); } + std::string KKMaterial::matElmDictionaryFileName() const { return findFile(elementsBaseName_); } + std::string KKMaterial::matIsoDictionaryFileName() const { return findFile(isotopesBaseName_); } + std::string KKMaterial::matMtrDictionaryFileName() const { return findFile(materialsBaseName_); } + } From 3477b3e767231717f953a96156a8bfb1a8e25473 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Sat, 16 May 2026 12:28:52 -0700 Subject: [PATCH 05/12] Move KKMaterial and KKStrawMaterial to KinKalGeom --- KinKalGeom/CMakeLists.txt | 4 ++++ {Mu2eKinKal => KinKalGeom}/inc/KKMaterial.hh | 8 ++++---- {Mu2eKinKal => KinKalGeom}/inc/KKStrawMaterial.hh | 0 {Mu2eKinKal => KinKalGeom}/src/KKMaterial.cc | 2 +- {Mu2eKinKal => KinKalGeom}/src/KKStrawMaterial.cc | 3 +-- KinKalGeom/src/SConscript | 3 +++ Mu2eKinKal/CMakeLists.txt | 2 -- Mu2eKinKal/inc/KKExtrap.hh | 2 +- Mu2eKinKal/inc/KKFit.hh | 2 +- Mu2eKinKal/inc/KKStrawXing.hh | 2 +- Mu2eKinKal/src/CentralHelixFit_module.cc | 2 +- Mu2eKinKal/src/KinematicLineFit_module.cc | 2 +- Mu2eKinKal/src/LoopHelixFit_module.cc | 2 +- Mu2eKinKal/src/RegrowKinematicLine_module.cc | 2 +- Mu2eKinKal/src/RegrowLoopHelix_module.cc | 2 +- RecoDataProducts/inc/TrkStraw.hh | 2 +- 16 files changed, 22 insertions(+), 18 deletions(-) rename {Mu2eKinKal => KinKalGeom}/inc/KKMaterial.hh (95%) rename {Mu2eKinKal => KinKalGeom}/inc/KKStrawMaterial.hh (100%) rename {Mu2eKinKal => KinKalGeom}/src/KKMaterial.cc (97%) rename {Mu2eKinKal => KinKalGeom}/src/KKStrawMaterial.cc (97%) diff --git a/KinKalGeom/CMakeLists.txt b/KinKalGeom/CMakeLists.txt index c600163619..7bc59e97ac 100644 --- a/KinKalGeom/CMakeLists.txt +++ b/KinKalGeom/CMakeLists.txt @@ -2,7 +2,11 @@ cet_make_library( SOURCE src/CRV.cc src/KinKalGeom.cc + src/KKMaterial.cc + src/KKStrawMaterial.cc LIBRARIES PUBLIC + Offline::ConfigTools + KinKal::MatEnv KinKal::Geometry KinKal::Trajectory KinKal::General diff --git a/Mu2eKinKal/inc/KKMaterial.hh b/KinKalGeom/inc/KKMaterial.hh similarity index 95% rename from Mu2eKinKal/inc/KKMaterial.hh rename to KinKalGeom/inc/KKMaterial.hh index fecefe8e04..f9b304b49f 100644 --- a/Mu2eKinKal/inc/KKMaterial.hh +++ b/KinKalGeom/inc/KKMaterial.hh @@ -1,5 +1,5 @@ -#ifndef Mu2eKinKal_KKMaterial_hh -#define Mu2eKinKal_KKMaterial_hh +#ifndef KinKalGeom_KKMaterial_hh +#define KinKalGeom_KKMaterial_hh // // build KinKal DetMaterial objects from art parameter configuration // @@ -9,8 +9,8 @@ // KinKal #include "KinKal/MatEnv/MatDBInfo.hh" #include "KinKal/MatEnv/FileFinderInterface.hh" -// Mu2eKinKal: moveme! -#include "Offline/Mu2eKinKal/inc/KKStrawMaterial.hh" +// KKGeom +#include "Offline/KinKalGeom/inc/KKStrawMaterial.hh" // mu2e #include "Offline/ConfigTools/inc/ConfigFileLookupPolicy.hh" diff --git a/Mu2eKinKal/inc/KKStrawMaterial.hh b/KinKalGeom/inc/KKStrawMaterial.hh similarity index 100% rename from Mu2eKinKal/inc/KKStrawMaterial.hh rename to KinKalGeom/inc/KKStrawMaterial.hh diff --git a/Mu2eKinKal/src/KKMaterial.cc b/KinKalGeom/src/KKMaterial.cc similarity index 97% rename from Mu2eKinKal/src/KKMaterial.cc rename to KinKalGeom/src/KKMaterial.cc index 72ff52ac6c..e80edefdb5 100644 --- a/Mu2eKinKal/src/KKMaterial.cc +++ b/KinKalGeom/src/KKMaterial.cc @@ -1,4 +1,4 @@ -#include "Offline/Mu2eKinKal/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/GeometryService/inc/GeomHandle.hh" #include "Offline/TrackerGeom/inc/Tracker.hh" #include "KinKal/MatEnv/DetMaterial.hh" diff --git a/Mu2eKinKal/src/KKStrawMaterial.cc b/KinKalGeom/src/KKStrawMaterial.cc similarity index 97% rename from Mu2eKinKal/src/KKStrawMaterial.cc rename to KinKalGeom/src/KKStrawMaterial.cc index 003d878fd3..22f0b808bb 100644 --- a/Mu2eKinKal/src/KKStrawMaterial.cc +++ b/KinKalGeom/src/KKStrawMaterial.cc @@ -1,5 +1,4 @@ -#include "Offline/Mu2eKinKal/inc/KKStrawMaterial.hh" -#include "Offline/Mu2eKinKal/inc/StrawXingUpdater.hh" +#include "Offline/KinKalGeom/inc/KKStrawMaterial.hh" #include "KinKal/MatEnv/DetMaterial.hh" #include "KinKal/MatEnv/MatDBInfo.hh" #include "KinKal/Detector/MaterialXing.hh" diff --git a/KinKalGeom/src/SConscript b/KinKalGeom/src/SConscript index decefd3fbe..325b8e8f5c 100644 --- a/KinKalGeom/src/SConscript +++ b/KinKalGeom/src/SConscript @@ -14,6 +14,8 @@ rootlibs = env['ROOTLIBS'] helper = mu2e_helper(env) mainlib = helper.make_mainlib([ + 'mu2e_ConfigTools', + 'KinKal_MatEnv', 'KinKal_Geometry', 'KinKal_Trajectory', 'KinKal_General', @@ -22,6 +24,7 @@ mainlib = helper.make_mainlib([ 'art_Framework_Services_Registry', 'art_Utilities', 'canvas', + 'cetlib', 'cetlib_except', 'CLHEP', 'Core', diff --git a/Mu2eKinKal/CMakeLists.txt b/Mu2eKinKal/CMakeLists.txt index 51fbad3f3e..3e1ae785bd 100644 --- a/Mu2eKinKal/CMakeLists.txt +++ b/Mu2eKinKal/CMakeLists.txt @@ -8,9 +8,7 @@ cet_make_library( src/KKConstantBField.cc src/KKFitSettings.cc src/KKFitUtilities.cc - src/KKMaterial.cc src/KKSHFlag.cc - src/KKStrawMaterial.cc src/StrawHitUpdaters.cc src/StrawXingUpdater.cc src/WHSIterator.cc diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index ab5a4205c7..6f67fc2a0c 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -13,7 +13,7 @@ #include "Offline/Mu2eKinKal/inc/ExtrapolateIPA.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateST.hh" #include "Offline/Mu2eKinKal/inc/KKShellXing.hh" -#include "Offline/Mu2eKinKal/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" #include "Offline/GeometryService/inc/GeomHandle.hh" #include "Offline/GeometryService/inc/GeometryService.hh" diff --git a/Mu2eKinKal/inc/KKFit.hh b/Mu2eKinKal/inc/KKFit.hh index 5b5077d35b..20c3a09b21 100644 --- a/Mu2eKinKal/inc/KKFit.hh +++ b/Mu2eKinKal/inc/KKFit.hh @@ -7,7 +7,7 @@ #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" -#include "Offline/Mu2eKinKal/inc/KKStrawMaterial.hh" +#include "Offline/KinKalGeom/inc/KKStrawMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKCaloHit.hh" #include "Offline/Mu2eKinKal/inc/KKFitUtilities.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" diff --git a/Mu2eKinKal/inc/KKStrawXing.hh b/Mu2eKinKal/inc/KKStrawXing.hh index e32f7e0239..1cc2136616 100644 --- a/Mu2eKinKal/inc/KKStrawXing.hh +++ b/Mu2eKinKal/inc/KKStrawXing.hh @@ -6,7 +6,7 @@ #include "KinKal/Detector/ElementXing.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/StrawXingUpdater.hh" -#include "Offline/Mu2eKinKal/inc/KKStrawMaterial.hh" +#include "Offline/KinKalGeom/inc/KKStrawMaterial.hh" #include "Offline/TrackerGeom/inc/Straw.hh" #include "Offline/DataProducts/inc/StrawId.hh" #include "KinKal/Trajectory/ParticleTrajectory.hh" diff --git a/Mu2eKinKal/src/CentralHelixFit_module.cc b/Mu2eKinKal/src/CentralHelixFit_module.cc index d45316dd45..eec7358e5d 100644 --- a/Mu2eKinKal/src/CentralHelixFit_module.cc +++ b/Mu2eKinKal/src/CentralHelixFit_module.cc @@ -52,7 +52,7 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" -#include "Offline/Mu2eKinKal/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" diff --git a/Mu2eKinKal/src/KinematicLineFit_module.cc b/Mu2eKinKal/src/KinematicLineFit_module.cc index 77cb6b665f..30e2079631 100644 --- a/Mu2eKinKal/src/KinematicLineFit_module.cc +++ b/Mu2eKinKal/src/KinematicLineFit_module.cc @@ -47,7 +47,7 @@ // Mu2eKinKal #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" -#include "Offline/Mu2eKinKal/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKBField.hh" #include "Offline/Mu2eKinKal/inc/KKFitUtilities.hh" diff --git a/Mu2eKinKal/src/LoopHelixFit_module.cc b/Mu2eKinKal/src/LoopHelixFit_module.cc index 1ac78588e8..d44d92fcb8 100644 --- a/Mu2eKinKal/src/LoopHelixFit_module.cc +++ b/Mu2eKinKal/src/LoopHelixFit_module.cc @@ -52,7 +52,7 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" -#include "Offline/Mu2eKinKal/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" diff --git a/Mu2eKinKal/src/RegrowKinematicLine_module.cc b/Mu2eKinKal/src/RegrowKinematicLine_module.cc index 451c6ecbe7..93d6e278d4 100644 --- a/Mu2eKinKal/src/RegrowKinematicLine_module.cc +++ b/Mu2eKinKal/src/RegrowKinematicLine_module.cc @@ -54,7 +54,7 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" -#include "Offline/Mu2eKinKal/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" diff --git a/Mu2eKinKal/src/RegrowLoopHelix_module.cc b/Mu2eKinKal/src/RegrowLoopHelix_module.cc index e4771e5646..2f648f2b6a 100644 --- a/Mu2eKinKal/src/RegrowLoopHelix_module.cc +++ b/Mu2eKinKal/src/RegrowLoopHelix_module.cc @@ -54,7 +54,7 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" -#include "Offline/Mu2eKinKal/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" diff --git a/RecoDataProducts/inc/TrkStraw.hh b/RecoDataProducts/inc/TrkStraw.hh index 59d4a6db99..306f2447e2 100644 --- a/RecoDataProducts/inc/TrkStraw.hh +++ b/RecoDataProducts/inc/TrkStraw.hh @@ -8,7 +8,7 @@ #define RecoDataProducts_TrkStraw_HH #include "Offline/DataProducts/inc/StrawId.hh" #include "Offline/RecoDataProducts/inc/StrawFlag.hh" -#include "Offline/Mu2eKinKal/inc/KKStrawMaterial.hh" +#include "Offline/KinKalGeom/inc/KKStrawMaterial.hh" #include "Offline/Mu2eKinKal/inc/StrawXingUpdater.hh" #include "KinKal/Trajectory/ClosestApproachData.hh" From e3dc0e0717d707b90d321b36da1d34e9f05b06bb Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Sat, 16 May 2026 14:05:17 -0700 Subject: [PATCH 06/12] Put KKMaterial into GeometryService --- GeometryService/inc/GeometryService.hh | 4 + GeometryService/inc/KinKalGeomMaker.hh | 1 + GeometryService/src/GeometryService.cc | 375 +++++++++++++------------ KinKalGeom/inc/KKMaterial.hh | 4 +- KinKalGeom/src/KKMaterial.cc | 2 +- Mu2eKinKal/src/SConscript | 4 +- fcl/standardServices.fcl | 14 + 7 files changed, 215 insertions(+), 189 deletions(-) diff --git a/GeometryService/inc/GeometryService.hh b/GeometryService/inc/GeometryService.hh index 806c80d2e2..f184ec1b02 100644 --- a/GeometryService/inc/GeometryService.hh +++ b/GeometryService/inc/GeometryService.hh @@ -21,6 +21,7 @@ #include "fhiclcpp/types/Table.h" #include "cetlib_except/exception.h" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/ConfigTools/inc/SimpleConfig.hh" #include "Offline/Mu2eInterfaces/inc/Detector.hh" #include "boost/shared_ptr.hpp" @@ -48,6 +49,7 @@ public: fhicl::Atom tool_type{Name("tool_type"),"Mu2e"}; }; + using KKMaterialConfig = KKMaterial::Config; struct Config { using Name=fhicl::Name; using Comment=fhicl::Comment; @@ -60,6 +62,7 @@ public: fhicl::Atom printConfig{Name("printConfig"),false}; fhicl::Atom printConfigTopLevel{Name("printConfigTopLevel"),false}; fhicl::Table simulatedDetector{Name("simulatedDetector")}; + fhicl::Table matSettings{Name("KinKalMaterial")}; }; using Parameters= art::ServiceTable; @@ -129,6 +132,7 @@ private: std::unique_ptr _bfConfig; const fhicl::ParameterSet _simulatedDetector; + const KKMaterialConfig _kkMat; // Load G4 geometry options std::unique_ptr _g4GeomOptions; diff --git a/GeometryService/inc/KinKalGeomMaker.hh b/GeometryService/inc/KinKalGeomMaker.hh index d0adbe8dbc..9c89168c9b 100644 --- a/GeometryService/inc/KinKalGeomMaker.hh +++ b/GeometryService/inc/KinKalGeomMaker.hh @@ -5,6 +5,7 @@ // Original author: Dave Brown (LBNL) 4/2026 // #include "Offline/KinKalGeom/inc/KinKalGeom.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" namespace mu2e { class KinKalGeomMaker { public: diff --git a/GeometryService/src/GeometryService.cc b/GeometryService/src/GeometryService.cc index a6c2786a90..91fe096b58 100644 --- a/GeometryService/src/GeometryService.cc +++ b/GeometryService/src/GeometryService.cc @@ -61,6 +61,7 @@ #include "Offline/BeamlineGeom/inc/TSdA.hh" #include "Offline/TrackerGeom/inc/Tracker.hh" #include "Offline/KinKalGeom/inc/KinKalGeom.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/CalorimeterGeom/inc/Calorimeter.hh" #include "Offline/GeometryService/inc/DiskCalorimeterMaker.hh" #include "Offline/CalorimeterGeom/inc/DiskCalorimeter.hh" @@ -96,7 +97,7 @@ using namespace std; namespace mu2e { GeometryService::GeometryService( const Parameters& pars, - art::ActivityRegistry&iRegistry) : + art::ActivityRegistry&iRegistry) : _inputfile( pars().inputFile()), _bFieldFile( pars().bFieldFile()), _allowReplacement( pars().allowReplacement()), @@ -107,6 +108,7 @@ namespace mu2e { _printTopLevel( pars().printConfigTopLevel()), _config(nullptr), _simulatedDetector( pars.get_PSet().get("simulatedDetector")), + _kkMat( pars().matSettings()), _standardMu2eDetector( _simulatedDetector.get("tool_type") == "Mu2e"), _detectors() { @@ -120,136 +122,136 @@ namespace mu2e { // This template can be defined here because this is a private method which is only // used by the code below in the same file. template - void GeometryService::addDetector(std::unique_ptr d) - { - if(_detectors.find(typeid(DET).name())!=_detectors.end()) { - throw cet::exception("GEOM") << "failed to install detector with type name " - << typeid(DET).name() << "\n"; - } + void GeometryService::addDetector(std::unique_ptr d) + { + if(_detectors.find(typeid(DET).name())!=_detectors.end()) { + throw cet::exception("GEOM") << "failed to install detector with type name " + << typeid(DET).name() << "\n"; + } DetectorPtr ptr(d.release()); _detectors[typeid(DET).name()] = ptr; - } + } template - void GeometryService::addDetectorAliasToBaseClass(std::unique_ptr d) - { + void GeometryService::addDetectorAliasToBaseClass(std::unique_ptr d) + { - std::string OriginalName = typeid(DET).name(); - DetMap::iterator it(_detectors.find(OriginalName)); + std::string OriginalName = typeid(DET).name(); + DetMap::iterator it(_detectors.find(OriginalName)); - if(it==_detectors.end()) - throw cet::exception("GEOM") - << "Can not alias an inexistant detector, detector " << OriginalName << "\n"; + if(it==_detectors.end()) + throw cet::exception("GEOM") + << "Can not alias an inexistant detector, detector " << OriginalName << "\n"; - std::string detectorName= typeid(DETALIAS).name() ; - _detectors[detectorName] = it->second; - } + std::string detectorName= typeid(DETALIAS).name() ; + _detectors[detectorName] = it->second; + } void - GeometryService::preBeginRun(art::Run const &) { - - if(++_run_count > 1) { - return; - } + GeometryService::preBeginRun(art::Run const &) { - _config = unique_ptr(new SimpleConfig(_inputfile, - _allowReplacement, - _messageOnReplacement, - _messageOnDefault )); - _config->printOpen(cout,"Geometry"); - - _bfConfig = unique_ptr(new SimpleConfig(_bFieldFile, - _allowReplacement, - _messageOnReplacement, - _messageOnDefault )); - _bfConfig->printOpen(cout,"BField"); - - - if(_printTopLevel) { - //print the top level geometry file contents - //the top level often contains a single named config file or a list of specific version files - ConfigFileLookupPolicy configFile; - std::string file = configFile(_inputfile); - std::ifstream in(file.c_str()); - if ( !in ) { - // No conf file for this test. - throw cet::exception("Geom") - << "GeometryService: Cannot open input file: " - << file - << endl; + if(++_run_count > 1) { + return; } - std::cout << "GeometryService: printing top level geometry file:\n"; - std::string line; - while ( in ){ - std::getline(in,line); - if ( !in ){ - break; - } - std::cout << line.c_str() << std::endl; + _config = unique_ptr(new SimpleConfig(_inputfile, + _allowReplacement, + _messageOnReplacement, + _messageOnDefault )); + _config->printOpen(cout,"Geometry"); + + _bfConfig = unique_ptr(new SimpleConfig(_bFieldFile, + _allowReplacement, + _messageOnReplacement, + _messageOnDefault )); + _bfConfig->printOpen(cout,"BField"); + + + if(_printTopLevel) { + //print the top level geometry file contents + //the top level often contains a single named config file or a list of specific version files + ConfigFileLookupPolicy configFile; + std::string file = configFile(_inputfile); + std::ifstream in(file.c_str()); + if ( !in ) { + // No conf file for this test. + throw cet::exception("Geom") + << "GeometryService: Cannot open input file: " + << file + << endl; + } + std::cout << "GeometryService: printing top level geometry file:\n"; + std::string line; + while ( in ){ + std::getline(in,line); + if ( !in ){ + break; + } + + std::cout << line.c_str() << std::endl; + } + std::cout << "GeometryService: finished printing top level geometry file.\n"; } - std::cout << "GeometryService: finished printing top level geometry file.\n"; - } - // Print final state of file after all substitutions. - if ( _printConfig ){ _config->print(cout, "Geom: "); } + // Print final state of file after all substitutions. + if ( _printConfig ){ _config->print(cout, "Geom: "); } - // 2019-03-24 P.M. : *not needed* decide if this is standard Mu2e detector or something else ... + // 2019-03-24 P.M. : *not needed* decide if this is standard Mu2e detector or something else ... - if (!isStandardMu2eDetector() || - !_config->getBool("mu2e.standardDetector",true)) { - cout << "Non-standard mu2e configuration, assuming it is intentional" << endl; - return; - } + if (!isStandardMu2eDetector() || + !_config->getBool("mu2e.standardDetector",true)) { + cout << "Non-standard mu2e configuration, assuming it is intentional" << endl; + return; + } - // Initialize geometry options - _g4GeomOptions = unique_ptr( new G4GeometryOptions( *_config ) ); + // Initialize geometry options + _g4GeomOptions = unique_ptr( new G4GeometryOptions( *_config ) ); - // Throw if the configuration is not self consistent. - checkConfig(); + // Throw if the configuration is not self consistent. + checkConfig(); - // This must be the first detector added since other makers may wish to use it. - std::unique_ptr tmpDetSys(DetectorSystemMaker::make(*_config)); - const DetectorSystem& detSys = *tmpDetSys.get(); - addDetector(std::move(tmpDetSys)); + // This must be the first detector added since other makers may wish to use it. + std::unique_ptr tmpDetSys(DetectorSystemMaker::make(*_config)); + const DetectorSystem& detSys = *tmpDetSys.get(); + addDetector(std::move(tmpDetSys)); - // Make a detector for every component present in the configuration. + // Make a detector for every component present in the configuration. - std::unique_ptr tmpBeamline(BeamlineMaker::make(*_config)); - const Beamline& beamline = *tmpBeamline.get(); - addDetector(std::move(tmpBeamline)); + std::unique_ptr tmpBeamline(BeamlineMaker::make(*_config)); + const Beamline& beamline = *tmpBeamline.get(); + addDetector(std::move(tmpBeamline)); - std::unique_ptr tmpProdTgt(ProductionTargetMaker::make(*_config, beamline.solenoidOffset())); - const ProductionTarget& prodTarget = *tmpProdTgt.get(); - addDetector(std::move(tmpProdTgt)); + std::unique_ptr tmpProdTgt(ProductionTargetMaker::make(*_config, beamline.solenoidOffset())); + const ProductionTarget& prodTarget = *tmpProdTgt.get(); + addDetector(std::move(tmpProdTgt)); - std::unique_ptr - tmpProductionSolenoid(ProductionSolenoidMaker(*_config, beamline.solenoidOffset()).getProductionSolenoidPtr()); + std::unique_ptr + tmpProductionSolenoid(ProductionSolenoidMaker(*_config, beamline.solenoidOffset()).getProductionSolenoidPtr()); - const ProductionSolenoid& ps = *tmpProductionSolenoid.get(); - addDetector(std::move(tmpProductionSolenoid)); + const ProductionSolenoid& ps = *tmpProductionSolenoid.get(); + addDetector(std::move(tmpProductionSolenoid)); - std::unique_ptr - tmpPSE(PSEnclosureMaker::make(*_config, ps.psEndRefPoint())); - const PSEnclosure& pse = *tmpPSE.get(); - addDetector(std::move(tmpPSE)); + std::unique_ptr + tmpPSE(PSEnclosureMaker::make(*_config, ps.psEndRefPoint())); + const PSEnclosure& pse = *tmpPSE.get(); + addDetector(std::move(tmpPSE)); - // The Z coordinate of the boundary between PS and TS vacua - StraightSection const * ts1vac = beamline.getTS().getTSVacuum( TransportSolenoid::TSRegion::TS1 ); - const double vacPS_TS_z = ts1vac->getGlobal().z() - ts1vac->getHalfLength(); + // The Z coordinate of the boundary between PS and TS vacua + StraightSection const * ts1vac = beamline.getTS().getTSVacuum( TransportSolenoid::TSRegion::TS1 ); + const double vacPS_TS_z = ts1vac->getGlobal().z() - ts1vac->getHalfLength(); - addDetector(PSVacuumMaker::make(*_config, ps, pse, vacPS_TS_z)); + addDetector(PSVacuumMaker::make(*_config, ps, pse, vacPS_TS_z)); - //addDetector(PSShieldMaker::make(*_config, ps.psEndRefPoint(), prodTarget.position())); + //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())); + 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 - 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;} @@ -257,118 +259,121 @@ namespace mu2e { - // Construct building solids - std::unique_ptr tmphall(Mu2eHallMaker::makeBuilding(*_g4GeomOptions,*_config)); - const Mu2eHall& hall = *tmphall.get(); + // Construct building solids + std::unique_ptr tmphall(Mu2eHallMaker::makeBuilding(*_g4GeomOptions,*_config)); + const Mu2eHall& hall = *tmphall.get(); - // Determine Mu2e envelope from building solids - std::unique_ptr mu2eEnv (new Mu2eEnvelope(hall,*_config)); + // Determine Mu2e envelope from building solids + std::unique_ptr mu2eEnv (new Mu2eEnvelope(hall,*_config)); - // Make dirt based on Mu2e envelope - Mu2eHallMaker::makeDirt( *tmphall.get(), *_g4GeomOptions, *_config, *mu2eEnv.get() ); - Mu2eHallMaker::makeRotated( *tmphall.get(), *_g4GeomOptions, *_config, *mu2eEnv.get() ); - Mu2eHallMaker::makeTrapDirt( *tmphall.get(), *_g4GeomOptions, *_config, *mu2eEnv.get() ); + // Make dirt based on Mu2e envelope + Mu2eHallMaker::makeDirt( *tmphall.get(), *_g4GeomOptions, *_config, *mu2eEnv.get() ); + Mu2eHallMaker::makeRotated( *tmphall.get(), *_g4GeomOptions, *_config, *mu2eEnv.get() ); + Mu2eHallMaker::makeTrapDirt( *tmphall.get(), *_g4GeomOptions, *_config, *mu2eEnv.get() ); - addDetector(std::move( tmphall ) ); - addDetector(std::move( mu2eEnv ) ); + addDetector(std::move( tmphall ) ); + addDetector(std::move( mu2eEnv ) ); - std::unique_ptr tmpDump(ProtonBeamDumpMaker::make(*_config, hall)); - const ProtonBeamDump& dump = *tmpDump.get(); - addDetector(std::move(tmpDump)); + std::unique_ptr tmpDump(ProtonBeamDumpMaker::make(*_config, hall)); + const ProtonBeamDump& dump = *tmpDump.get(); + addDetector(std::move(tmpDump)); - // beamline info used to position DS - std::unique_ptr tmpDS( DetectorSolenoidMaker::make( *_config, beamline ) ); - const DetectorSolenoid& ds = *tmpDS.get(); - addDetector(std::move(tmpDS)); + // beamline info used to position DS + std::unique_ptr tmpDS( DetectorSolenoidMaker::make( *_config, beamline ) ); + const DetectorSolenoid& ds = *tmpDS.get(); + addDetector(std::move(tmpDS)); - // DS info used to position DS downstream shielding - addDetector( DetectorSolenoidShieldingMaker::make( *_config, ds ) ); + // DS info used to position DS downstream shielding + addDetector( DetectorSolenoidShieldingMaker::make( *_config, ds ) ); - std::unique_ptr tmptgt(StoppingTargetMaker(detSys.getOrigin(), *_config).getTargetPtr()); - const StoppingTarget& target = *tmptgt.get(); - addDetector(std::move(tmptgt)); + std::unique_ptr tmptgt(StoppingTargetMaker(detSys.getOrigin(), *_config).getTargetPtr()); + const StoppingTarget& target = *tmptgt.get(); + addDetector(std::move(tmptgt)); - if (_config->getBool("hasTracker",false)){ - TrackerMaker ttm( *_config ); - addDetector( ttm.getTrackerPtr() ); - } + if (_config->getBool("hasTracker",false)){ + TrackerMaker ttm( *_config ); + addDetector( ttm.getTrackerPtr() ); + } - if(_config->getBool("hasMBS",false)){ - MBSMaker mbs( *_config, beamline.solenoidOffset() ); - addDetector( mbs.getMBSPtr() ); - } + if(_config->getBool("hasMBS",false)){ + MBSMaker mbs( *_config, beamline.solenoidOffset() ); + addDetector( mbs.getMBSPtr() ); + } - if(_config->getBool("hasDiskCalorimeter",false)){ - DiskCalorimeterMaker calorm( *_config, beamline.solenoidOffset() ); - addDetector( calorm.calorimeterPtr() ); - addDetectorAliasToBaseClass( calorm.calorimeterPtr() ); //add an alias to detector list - } + if(_config->getBool("hasDiskCalorimeter",false)){ + DiskCalorimeterMaker calorm( *_config, beamline.solenoidOffset() ); + addDetector( calorm.calorimeterPtr() ); + addDetectorAliasToBaseClass( calorm.calorimeterPtr() ); //add an alias to detector list + } - if(_config->getBool("hasCosmicRayShield",false)){ - CosmicRayShieldMaker crs( *_config, beamline.solenoidOffset() ); - addDetector( crs.getCosmicRayShieldPtr() ); - } + if(_config->getBool("hasCosmicRayShield",false)){ + CosmicRayShieldMaker crs( *_config, beamline.solenoidOffset() ); + addDetector( crs.getCosmicRayShieldPtr() ); + } - if(_config->getBool("hasTSdA",false)){ - addDetector( TSdAMaker::make(*_config,ds) ); - } + if(_config->getBool("hasTSdA",false)){ + addDetector( TSdAMaker::make(*_config,ds) ); + } - if(_config->getBool("hasExternalShielding",false)) { - addDetector( ExtShieldUpstreamMaker::make(*_config) ); - addDetector( ExtShieldDownstreamMaker::make(*_config)); - addDetector( SaddleMaker::make(*_config)); - addDetector( PipeMaker::make(*_config)); - addDetector( ElectronicRackMaker::make(*_config)); - } + if(_config->getBool("hasExternalShielding",false)) { + addDetector( ExtShieldUpstreamMaker::make(*_config) ); + addDetector( ExtShieldDownstreamMaker::make(*_config)); + addDetector( SaddleMaker::make(*_config)); + addDetector( PipeMaker::make(*_config)); + addDetector( ElectronicRackMaker::make(*_config)); + } - std::unique_ptr tmpemb(ExtMonFNALBuildingMaker::make(*_config, hall, dump)); - const ExtMonFNALBuilding& emfb = *tmpemb.get(); - addDetector(std::move(tmpemb)); - if(_config->getBool("hasExtMonFNAL",false)){ - addDetector(ExtMonFNAL::ExtMonMaker::make(*_config, emfb)); - } + std::unique_ptr tmpemb(ExtMonFNALBuildingMaker::make(*_config, hall, dump)); + const ExtMonFNALBuilding& emfb = *tmpemb.get(); + addDetector(std::move(tmpemb)); + if(_config->getBool("hasExtMonFNAL",false)){ + addDetector(ExtMonFNAL::ExtMonMaker::make(*_config, emfb)); + } - if (_config->getBool("hasPTM",false) ){ - std::unique_ptr ptmon(PTMMaker::make(*_config)); - addDetector(std::move(ptmon)); - } + if (_config->getBool("hasPTM",false) ){ + std::unique_ptr ptmon(PTMMaker::make(*_config)); + addDetector(std::move(ptmon)); + } - if(_config->getBool("hasSTM",false)){ - STMMaker stm( *_config, beamline.solenoidOffset() ); - addDetector( stm.getSTMPtr() ); - } + if(_config->getBool("hasSTM",false)){ + STMMaker stm( *_config, beamline.solenoidOffset() ); + addDetector( stm.getSTMPtr() ); + } - if(_config->getBool("hasVirtualDetector",false)){ - addDetector(VirtualDetectorMaker::make(*_config)); - } + if(_config->getBool("hasVirtualDetector",false)){ + addDetector(VirtualDetectorMaker::make(*_config)); + } - if(_bfConfig->getBool("hasBFieldManager",false)){ - std::unique_ptr bfc( BFieldConfigMaker(*_bfConfig, beamline).getBFieldConfig() ); - BFieldManagerMaker bfmgr(*bfc); - addDetector(std::move(bfc)); - addDetector(bfmgr.getBFieldManager()); - } + if(_bfConfig->getBool("hasBFieldManager",false)){ + std::unique_ptr bfc( BFieldConfigMaker(*_bfConfig, beamline).getBFieldConfig() ); + BFieldManagerMaker bfmgr(*bfc); + addDetector(std::move(bfc)); + addDetector(bfmgr.getBFieldManager()); + } - if(_config->getBool("hasProtonAbsorber",false) && !_config->getBool("protonabsorber.isHelical", false) ){ - MECOStyleProtonAbsorberMaker mecopam( *_config, ds, target); - addDetector( mecopam.getMECOStyleProtonAbsorberPtr() ); - } + if(_config->getBool("hasProtonAbsorber",false) && !_config->getBool("protonabsorber.isHelical", false) ){ + MECOStyleProtonAbsorberMaker mecopam( *_config, ds, target); + addDetector( mecopam.getMECOStyleProtonAbsorberPtr() ); + } - // This class has a default c'tor with all available information internally. - std::unique_ptr dusafMu2e{ std::make_unique() }; - addDetector( std::move(dusafMu2e) ); + // This class has a default c'tor with all available information internally. + std::unique_ptr dusafMu2e{ std::make_unique() }; + addDetector( std::move(dusafMu2e) ); - // build KinKalGeom, used in track reconstruction and extrapolation - KinKalGeomMaker kkgm; - addDetector( std::move(kkgm.makeKKG()) ); + // build KinKalGeom, used in track reconstruction and extrapolation + KinKalGeomMaker kkgm; + addDetector( std::move(kkgm.makeKKG()) ); + // directly build KKMaterial; it's constructor does everything + addDetector( std::make_unique(_kkMat)); + return ; - } // preBeginRun() + } // preBeginRun() // Check that the configuration is self consistent. void GeometryService::checkConfig(){ diff --git a/KinKalGeom/inc/KKMaterial.hh b/KinKalGeom/inc/KKMaterial.hh index f9b304b49f..13b83d62eb 100644 --- a/KinKalGeom/inc/KKMaterial.hh +++ b/KinKalGeom/inc/KKMaterial.hh @@ -13,12 +13,14 @@ #include "Offline/KinKalGeom/inc/KKStrawMaterial.hh" // mu2e #include "Offline/ConfigTools/inc/ConfigFileLookupPolicy.hh" +#include "Offline/Mu2eInterfaces/inc/Detector.hh" +#include "Offline/Mu2eInterfaces/inc/ProditionsEntity.hh" #include #include namespace mu2e { - class KKMaterial : public MatEnv::FileFinderInterface { + class KKMaterial : public MatEnv::FileFinderInterface, public Detector, public ProditionsEntity { public: using Name = fhicl::Name; using Comment = fhicl::Comment; diff --git a/KinKalGeom/src/KKMaterial.cc b/KinKalGeom/src/KKMaterial.cc index e80edefdb5..8c390445d6 100644 --- a/KinKalGeom/src/KKMaterial.cc +++ b/KinKalGeom/src/KKMaterial.cc @@ -7,7 +7,7 @@ namespace mu2e { using MatDBInfo = MatEnv::MatDBInfo; using MatEnv::DetMaterial; - KKMaterial::KKMaterial(KKMaterial::Config const& matconfig) : + KKMaterial::KKMaterial(KKMaterial::Config const& matconfig) : ProditionsEntity("KKMaterial"), elementsBaseName_(matconfig.elements()), isotopesBaseName_(matconfig.isotopes()), materialsBaseName_(matconfig.materials()), diff --git a/Mu2eKinKal/src/SConscript b/Mu2eKinKal/src/SConscript index 33a145e8e5..1419c781cc 100644 --- a/Mu2eKinKal/src/SConscript +++ b/Mu2eKinKal/src/SConscript @@ -20,6 +20,7 @@ mainlib = helper.make_mainlib([ 'mu2e_GeometryService', 'mu2e_BFieldGeom', 'mu2e_CalorimeterGeom', + 'mu2e_KinKalGeom', 'mu2e_TrackerGeom', 'mu2e_RecoDataProducts', 'mu2e_GlobalConstantsService', @@ -40,7 +41,6 @@ mainlib = helper.make_mainlib([ 'tbb', 'cetlib', 'cetlib_except', - 'mu2e_KinKalGeom', 'KinKal_Fit', 'KinKal_Trajectory', 'KinKal_Geometry', @@ -62,6 +62,7 @@ helper.make_plugins([mainlib, 'mu2e_GeometryService', 'mu2e_BFieldGeom', 'mu2e_CalorimeterGeom', + 'mu2e_KinKalGeom', 'mu2e_TrackerGeom', 'mu2e_RecoDataProducts', 'mu2e_GlobalConstantsService', @@ -87,7 +88,6 @@ helper.make_plugins([mainlib, 'cetlib', 'cetlib_except', 'CLHEP', - 'mu2e_KinKalGeom', 'KinKal_Fit', 'KinKal_Trajectory', 'KinKal_Geometry', diff --git a/fcl/standardServices.fcl b/fcl/standardServices.fcl index bbbc18dbab..96191e747b 100644 --- a/fcl/standardServices.fcl +++ b/fcl/standardServices.fcl @@ -42,6 +42,20 @@ Services : { inputFile : "Offline/Mu2eG4/geom/geom_common.txt" bFieldFile : "Offline/Mu2eG4/geom/bfgeom_v01.txt" simulatedDetector : { tool_type: "Mu2e" } + KinKalMaterial : { + elements : "Offline/TrackerConditions/data/ElementsList.data" + isotopes : "Offline/TrackerConditions/data/IsotopesList.data" + materials : "Offline/TrackerConditions/data/MaterialsList.data" + strawGasMaterialName : "straw-gas" + strawWallMaterialName : "straw-wall" + strawWireMaterialName : "straw-wire" + IPAMaterialName : "HDPE" + STMaterialName : "Target" + IonizationEnergyLossMode : 1 # Moyal mean + SolidScatteringFraction : 0.999999 + GasScatteringFraction : 0.9999999 + ElectronBrehmsFraction : 0.04 + } } GlobalConstantsService : { inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt" } DbService : @local::DbEmpty From 3dbe69724abb1b97a0d8a7b7dc94ca2546574cc9 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Sat, 16 May 2026 14:35:20 -0700 Subject: [PATCH 07/12] Convert to using KKMaterial from GeometryService --- Mu2eKinKal/fcl/prolog.fcl | 20 -------- Mu2eKinKal/inc/KKExtrap.hh | 49 +++++++++----------- Mu2eKinKal/src/CentralHelixFit_module.cc | 9 ++-- Mu2eKinKal/src/KinematicLineFit_module.cc | 17 +++---- Mu2eKinKal/src/LoopHelixFit_module.cc | 15 ++---- Mu2eKinKal/src/RegrowKinematicLine_module.cc | 12 ++--- Mu2eKinKal/src/RegrowLoopHelix_module.cc | 10 ++-- 7 files changed, 43 insertions(+), 89 deletions(-) diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index c827389d0f..808ff07c56 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -3,20 +3,6 @@ BEGIN_PROLOG Mu2eKinKal : { # general configuration # - MAT: { - elements : "Offline/TrackerConditions/data/ElementsList.data" - isotopes : "Offline/TrackerConditions/data/IsotopesList.data" - materials : "Offline/TrackerConditions/data/MaterialsList.data" - strawGasMaterialName : "straw-gas" - strawWallMaterialName : "straw-wall" - strawWireMaterialName : "straw-wire" - IPAMaterialName : "HDPE" - STMaterialName : "Target" - IonizationEnergyLossMode : 1 # Moyal mean - SolidScatteringFraction : 0.999999 - GasScatteringFraction : 0.9999999 - ElectronBrehmsFraction : 0.04 - } KKFIT: { PrintLevel : 0 TPOCAPrecision : 1e-4 # mm @@ -446,7 +432,6 @@ Mu2eKinKal : { LHSeedFit : { module_type : LoopHelixFit - MaterialSettings : @local::Mu2eKinKal.MAT KKFitSettings: { @table::Mu2eKinKal.KKFIT SampleSurfaces : ["TT_Mid"] @@ -468,7 +453,6 @@ Mu2eKinKal : { LHDriftFit : { module_type : LoopHelixFit - MaterialSettings : @local::Mu2eKinKal.MAT KKFitSettings : { @table::Mu2eKinKal.KKFIT SampleSurfaces : ["ST_Outer","ST_Front","ST_Back"] # these are additional surfaces; surfaces used in extrapolation are also sampled @@ -493,7 +477,6 @@ Mu2eKinKal : { KLSeedFit : { module_type : KinematicLineFit - MaterialSettings : @local::Mu2eKinKal.MAT KKFitSettings: @local::Mu2eKinKal.KKFIT FitSettings : @local::Mu2eKinKal.CHSEEDFIT ExtensionSettings : { @@ -508,7 +491,6 @@ Mu2eKinKal : { KLDriftFit : { module_type : KinematicLineFit - MaterialSettings : @local::Mu2eKinKal.MAT KKFitSettings: { @table::Mu2eKinKal.KKFIT # DNB I don't know where this time offset difference WRT helical fits comes from. If it's physical, we will need a better @@ -531,7 +513,6 @@ Mu2eKinKal : { CHSeedFit : { module_type : CentralHelixFit - MaterialSettings : @local::Mu2eKinKal.MAT KKFitSettings: { @table::Mu2eKinKal.KKFIT SaveTrajectory : T0 @@ -547,7 +528,6 @@ Mu2eKinKal : { CHDriftFit : { module_type : CentralHelixFit - MaterialSettings : @local::Mu2eKinKal.MAT KKFitSettings: { @table::Mu2eKinKal.KKFIT SaveTrajectory : Detector diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index 6f67fc2a0c..f46290275d 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -30,7 +30,7 @@ namespace mu2e { class KKExtrap { public: - explicit KKExtrap(KKExtrapConfig const& exconfig,KKMaterial const& kkmat); + explicit KKExtrap(KKExtrapConfig const& exconfig); // extrapolation functions; these are templated on the type of trajectory template void extrapolate(KKTrack& ktrk) const; template void toTrackerEnds(KKTrack& ktrk) const; @@ -43,18 +43,16 @@ namespace mu2e { private: int debug_; double btol_, intertol_, maxdt_; - KKMaterial const& kkmat_; bool backToTracker_, toOPA_, toTrackerEnds_, upstream_; double ipathick_ = 0.511; // ipa thickness: should come from geometry service TODO double stthick_ = 0.1056; // st foil thickness: should come from geometry service TODO }; - KKExtrap::KKExtrap(KKExtrapConfig const& extrapconfig,KKMaterial const& kkmat) : + KKExtrap::KKExtrap(KKExtrapConfig const& extrapconfig) : debug_(extrapconfig.Debug()), btol_(extrapconfig.btol()), intertol_(extrapconfig.interTol()), maxdt_(extrapconfig.MaxDt()), - kkmat_(kkmat), backToTracker_(extrapconfig.BackToTracker()), toOPA_(extrapconfig.ToOPA()), toTrackerEnds_(extrapconfig.ToTrackerEnds()), @@ -63,8 +61,7 @@ namespace mu2e { template void KKExtrap::extrapolate(KKTrack& ktrk) const { GeomHandle kkg_h; - auto const& kkg = *kkg_h; - ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker()->front().center().Z(),debug_); + ExtrapolateToZ trackerFront(maxdt_,btol_,kkg_h->tracker()->front().center().Z(),debug_); // define the time direction according to the fit direction inside the tracker auto const& ftraj = ktrk.fitTraj(); @@ -103,10 +100,9 @@ namespace mu2e { template void KKExtrap::toTrackerEnds(KKTrack& ktrk) const { GeomHandle kkg_h; - auto const& kkg = *kkg_h; - ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker()->front().center().Z(),debug_); - ExtrapolateToZ trackerBack(maxdt_,btol_,kkg.tracker()->back().center().Z(),debug_); + ExtrapolateToZ trackerFront(maxdt_,btol_,kkg_h->tracker()->front().center().Z(),debug_); + ExtrapolateToZ trackerBack(maxdt_,btol_,kkg_h->tracker()->back().center().Z(),debug_); // time direction to reach the bounding surfaces from the active region depends on the z momentum. This calculation assumes the particle doesn't // reflect inside the tracker volume auto const& ftraj = ktrk.fitTraj(); @@ -121,18 +117,18 @@ namespace mu2e { static const SurfaceId tt_back("TT_Back"); // start with the middle - auto midinter = KinKal::intersect(ftraj,kkg.tracker()->middle(),ftraj.range(),intertol_); + auto midinter = KinKal::intersect(ftraj,kkg_h->tracker()->middle(),ftraj.range(),intertol_); if(midinter.good()) ktrk.addIntersection(tt_mid,midinter); if(tofront){ // check the front piece first; that is usually correct // track extrapolation to the front succeeded, but the intersection failed. Use the last trajectory to force an intersection auto fhel = fronttdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); - auto frontinter = KinKal::intersect(fhel,kkg.tracker()->front(),fhel.range(),intertol_,fronttdir); + auto frontinter = KinKal::intersect(fhel,kkg_h->tracker()->front(),fhel.range(),intertol_,fronttdir); if(!frontinter.good()){ // start from the middle TimeRange frange = ftraj.range(); if(midinter.good())frange = fronttdir == TimeDir::forwards ? TimeRange(midinter.time_,ftraj.range().end()) : TimeRange(ftraj.range().begin(),midinter.time_); - frontinter = KinKal::intersect(ftraj,kkg.tracker()->front(),frange,intertol_,fronttdir); + frontinter = KinKal::intersect(ftraj,kkg_h->tracker()->front(),frange,intertol_,fronttdir); } if(frontinter.good()) ktrk.addIntersection(tt_front,frontinter); } @@ -140,19 +136,19 @@ namespace mu2e { // start from the middle TimeRange brange = ftraj.range(); if(midinter.good())brange = backtdir == TimeDir::forwards ? TimeRange(midinter.time_,ftraj.range().end()) : TimeRange(ftraj.range().begin(),midinter.time_); - auto backinter = KinKal::intersect(ftraj,kkg.tracker()->back(),brange,intertol_,backtdir); + auto backinter = KinKal::intersect(ftraj,kkg_h->tracker()->back(),brange,intertol_,backtdir); if(backinter.good())ktrk.addIntersection(tt_back,backinter); } } template bool KKExtrap::extrapolateIPA(KKTrack& ktrk,TimeDir tdir) const { GeomHandle kkg_h; - auto const& kkg = *kkg_h; + GeomHandle kkmat_h; using KKIPAXING = KKShellXing; using KKIPAXINGPTR = std::shared_ptr; // extraplate the fit through the IPA. This will add material effects for each intersection. It will continue till the // track exits the IPA - ExtrapolateIPA extrapIPA(maxdt_,btol_,intertol_,kkg.DS()->innerProtonAbsorberPtr(),debug_); + ExtrapolateIPA extrapIPA(maxdt_,btol_,intertol_,kkg_h->DS()->innerProtonAbsorberPtr(),debug_); if(extrapIPA.debug() > 0)std::cout << "extrapolating to IPA " << std::endl; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId IPASID("IPA"); @@ -163,8 +159,8 @@ namespace mu2e { if(extrapIPA.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(); - auto const& IPA = kkg.DS()->innerProtonAbsorberPtr(); - KKIPAXINGPTR ipaxingptr = std::make_shared(IPA,IPASID,*kkmat_.IPAMaterial(),extrapIPA.intersection(),reftrajptr,ipathick_,extrapIPA.interTolerance()); + auto const& IPA = kkg_h->DS()->innerProtonAbsorberPtr(); + KKIPAXINGPTR ipaxingptr = std::make_shared(IPA,IPASID,*kkmat_h->IPAMaterial(),extrapIPA.intersection(),reftrajptr,ipathick_,extrapIPA.interTolerance()); if(extrapIPA.debug() > 0){ double dmom, paramomvar, perpmomvar; ipaxingptr->materialEffects(dmom,paramomvar,perpmomvar); @@ -191,11 +187,11 @@ namespace mu2e { using KKSTXING = KKShellXing; using KKSTXINGPTR = std::shared_ptr; GeomHandle kkg_h; - auto const& kkg = *kkg_h; + GeomHandle kkmat_h; // extraplate the fit through the ST. This will add material effects for each foil intersection. It will continue till the // track exits the ST in Z - ExtrapolateST extrapST(maxdt_,btol_,intertol_,*kkg.ST(),debug_); + ExtrapolateST extrapST(maxdt_,btol_,intertol_,*kkg_h->ST(),debug_); auto const& ftraj = ktrk.fitTraj(); double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); auto startdir = ftraj.direction(starttime); @@ -205,7 +201,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_.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()); if(extrapST.debug() > 0){ double dmom, paramomvar, perpmomvar; stxingptr->materialEffects(dmom,paramomvar,perpmomvar); @@ -231,15 +227,14 @@ namespace mu2e { template bool KKExtrap::extrapolateTracker(KKTrack& ktrk,TimeDir tdir) const { GeomHandle kkg_h; - auto const& kkg = *kkg_h; - ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker()->front().center().Z(),debug_); + ExtrapolateToZ trackerFront(maxdt_,btol_,kkg_h->tracker()->front().center().Z(),debug_); if(trackerFront.debug() > 0)std::cout << "extrapolating to Tracker " << std::endl; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId TrackerSID("TT_Front"); ktrk.extrapolate(tdir,trackerFront); // the last piece appended should cover the necessary range auto const& ktraj = tdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); - auto trkfrontinter = KinKal::intersect(ftraj,kkg.tracker()->front(),ktraj.range(),intertol_,tdir); + auto trkfrontinter = KinKal::intersect(ftraj,kkg_h->tracker()->front(),ktraj.range(),intertol_,tdir); if(trkfrontinter.onsurface_){ // dont worry about bounds here ktrk.addIntersection(TrackerSID,trkfrontinter); return true; @@ -249,8 +244,7 @@ namespace mu2e { template bool KKExtrap::extrapolateTSDA(KKTrack& ktrk,TimeDir tdir) const { GeomHandle kkg_h; - auto const& kkg = *kkg_h; - ExtrapolateToZ TSDA(maxdt_,btol_,kkg.DS()->upstreamAbsorber().center().Z(),debug_); + ExtrapolateToZ TSDA(maxdt_,btol_,kkg_h->DS()->upstreamAbsorber().center().Z(),debug_); if(TSDA.debug() > 0)std::cout << "extrapolating to TSDA " << std::endl; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId TSDASID("TSDA"); @@ -261,7 +255,7 @@ namespace mu2e { bool retval = epos.Z() < TSDA.zVal(); if(retval){ auto const& ktraj = tdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); - auto tsdainter = KinKal::intersect(ftraj,kkg.DS()->upstreamAbsorber(),ktraj.range(),intertol_,tdir); + auto tsdainter = KinKal::intersect(ftraj,kkg_h->DS()->upstreamAbsorber(),ktraj.range(),intertol_,tdir); if(tsdainter.onsurface_)ktrk.addIntersection(TSDASID,tsdainter); } return retval; @@ -269,11 +263,10 @@ namespace mu2e { template void KKExtrap::toOPA(KKTrack& ktrk, double tstart, TimeDir tdir) const { GeomHandle kkg_h; - auto const& kkg = *kkg_h; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId OPASID("OPA"); TimeRange trange = tdir == TimeDir::forwards ? TimeRange(tstart,ftraj.range().end()) : TimeRange(ftraj.range().begin(),tstart); - auto opainter = KinKal::intersect(ftraj,kkg.DS()->outerProtonAbsorber(),trange,intertol_,tdir); + auto opainter = KinKal::intersect(ftraj,kkg_h->DS()->outerProtonAbsorber(),trange,intertol_,tdir); if(opainter.good()){ ktrk.addIntersection(OPASID,opainter); } diff --git a/Mu2eKinKal/src/CentralHelixFit_module.cc b/Mu2eKinKal/src/CentralHelixFit_module.cc index eec7358e5d..f6cef066ed 100644 --- a/Mu2eKinKal/src/CentralHelixFit_module.cc +++ b/Mu2eKinKal/src/CentralHelixFit_module.cc @@ -102,7 +102,6 @@ namespace mu2e { using KKFitConfig = Mu2eKinKal::KKFitConfig; using KKModuleConfig = Mu2eKinKal::KKModuleConfig; - using KKMaterialConfig = KKMaterial::Config; using Name = fhicl::Name; using Comment = fhicl::Comment; using KinKal::DVEC; @@ -128,7 +127,6 @@ namespace mu2e { fhicl::Table kkfitSettings { Name("KKFitSettings") }; fhicl::Table fitSettings { Name("FitSettings") }; fhicl::Table extSettings { Name("ExtensionSettings") }; - fhicl::Table matSettings { Name("MaterialSettings") }; // helix module specific config }; @@ -155,7 +153,6 @@ namespace mu2e { int print_; PDGCode::type fpart_; KKFIT kkfit_; // fit helper - KKMaterial kkmat_; // material helper DMAT seedcov_; // seed covariance matrix double mass_; // particle mass int PDGcharge_; // PDG particle charge @@ -184,7 +181,6 @@ namespace mu2e { print_(settings().modSettings().printLevel()), fpart_(static_cast(settings().modSettings().fitParticle())), kkfit_(settings().kkfitSettings()), - kkmat_(settings().matSettings()), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), exconfig_(Mu2eKinKal::makeConfig(settings().extSettings())), fixedfield_(false), @@ -246,6 +242,7 @@ namespace mu2e { void CentralHelixFit::produce(art::Event& event ) { GeomHandle calo_h; GeomHandle nominalTracker_h; + GeomHandle kkmat_h; // find current proditions auto const& strawresponse = strawResponse_h_.getPtr(event.id()); auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); @@ -319,7 +316,7 @@ namespace mu2e { strawhits.reserve(strawHitIdxs.size()); KKSTRAWXINGCOL strawxings; strawxings.reserve(strawHitIdxs.size()); - kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_.strawMaterial(), pseedtraj, chcol, strawHitIdxs, strawhits, strawxings); + kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_h->strawMaterial(), pseedtraj, chcol, strawHitIdxs, strawhits, strawxings); // optionally (and if present) add the CaloCluster as a constraint // verify the cluster looks physically reasonable before adding it TODO! Or, let the KKCaloHit updater do it TODO KKCALOHITCOL calohits; @@ -335,7 +332,7 @@ namespace mu2e { // if we have an extension schedule, extend. if(goodfit && exconfig_.schedule().size() > 0) { // std::cout << "EXTENDING TRACK " << event.id() << " " << index << std::endl; - kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_.strawMaterial(), chcol, *calo_h, cc_H, *kktrk ); + kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H, *kktrk ); goodfit = goodFit(*kktrk); } diff --git a/Mu2eKinKal/src/KinematicLineFit_module.cc b/Mu2eKinKal/src/KinematicLineFit_module.cc index 30e2079631..60f135404c 100644 --- a/Mu2eKinKal/src/KinematicLineFit_module.cc +++ b/Mu2eKinKal/src/KinematicLineFit_module.cc @@ -103,7 +103,6 @@ namespace mu2e { using KKConfig = Mu2eKinKal::KinKalConfig; using KKFitConfig = Mu2eKinKal::KKFitConfig; using KKModuleConfig = Mu2eKinKal::KKModuleConfig; - using KKMaterialConfig = KKMaterial::Config; class KinematicLineFit : public art::EDProducer { using Name = fhicl::Name; @@ -133,7 +132,6 @@ namespace mu2e { fhicl::Table mu2eSettings { Name("KKFitSettings") }; fhicl::Table fitSettings { Name("FitSettings") }; fhicl::Table extSettings { Name("ExtensionSettings") }; - fhicl::Table matSettings { Name("MaterialSettings") }; fhicl::OptionalTable Extrapolation { Name("Extrapolation") }; fhicl::OptionalAtom fitDirection { Name("FitDirection"), Comment("Particle direction to fit, either \"upstream\" or \"downstream\"")}; }; @@ -163,7 +161,6 @@ namespace mu2e { PDGCode::type fpart_; TrkFitDirection fdir_; KKFIT kkfit_; // fit helper - KKMaterial kkmat_; // material helper DMAT seedcov_; // seed covariance matrix std::array paramconstraints_; double mass_; // particle mass @@ -191,7 +188,6 @@ namespace mu2e { seedmom_(settings().modSettings().seedmom()), fpart_(static_cast(settings().modSettings().fitParticle())), kkfit_(settings().mu2eSettings()), - kkmat_(settings().matSettings()), intertol_(settings().modSettings().interTol()), sampletbuff_(settings().modSettings().sampleTBuff()), sampleinrange_(settings().modSettings().sampleInRange()), @@ -257,8 +253,7 @@ namespace mu2e { kkbf_ = std::make_unique(*bfmgr,*det); // translate the sample surface names to actual surfaces using the KinKalGeom. This must be done after construction as the KKGeom object now comes from GeometryService GeomHandle kkg_h; - auto const& kkg = *kkg_h; - kkg.surfaces(ssids_,surfacess_to_sample_); + kkg_h->surfaces(ssids_,surfacess_to_sample_); } void KinematicLineFit::produce(art::Event& event ) { @@ -303,7 +298,7 @@ namespace mu2e { KKSTRAWXINGCOL strawxings; strawhits.reserve(strawHitIdxs.size()); strawxings.reserve(strawHitIdxs.size()); - kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_.strawMaterial(), pseedtraj, *chcolptr, strawHitIdxs, strawhits, strawxings); + kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_h->strawMaterial(), pseedtraj, *chcolptr, strawHitIdxs, strawhits, strawxings); //here KKCALOHITCOL calohits; @@ -324,7 +319,7 @@ namespace mu2e { auto kktrk = make_unique(config_,*kkbf_,seedtraj,fpart_,kkfit_.strawHitClusterer(),strawhits,strawxings,calohits,paramconstraints_); auto goodfit = goodFit(*kktrk); if(goodfit && exconfig_.schedule().size() > 0){ - kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_.strawMaterial(), chcol, *calo_h, cc_H, *kktrk ); + kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H, *kktrk ); } goodfit = goodFit(*kktrk); // extrapolate as required @@ -412,9 +407,9 @@ namespace mu2e { void KinematicLineFit::extrapolate(KKTRK& ktrk) const { GeomHandle kkg_h; - auto const& kkg = *kkg_h; + GeomHandle kkmat_h; // extrapolate to the extracted CRV. This function should be migrated to KKExtrap TODO - auto TCRV = ExtrapolateTCRV(maxdt_,btol_,intertol_,minv_,*kkg.TCRV(),extrapdebug_); + auto TCRV = ExtrapolateTCRV(maxdt_,btol_,intertol_,minv_,*kkg_h->TCRV(),extrapdebug_); auto const& ftraj = ktrk.fitTraj(); static const SurfaceId TCRVSID("TCRV"); @@ -437,7 +432,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(); // TODO add DS and shielding material - KKCRVXINGPTR crvxingptr = std::make_shared(TCRV.module(), TCRVSID, *kkmat_.STMaterial(),TCRV.intersection(),reftrajptr,tcrvthick_,TCRV.interTolerance()); + KKCRVXINGPTR crvxingptr = std::make_shared(TCRV.module(), TCRVSID, *kkmat_h->STMaterial(),TCRV.intersection(),reftrajptr,tcrvthick_,TCRV.interTolerance()); ktrk.addTCRVXing(crvxingptr,tdir); } } while(hadintersection); diff --git a/Mu2eKinKal/src/LoopHelixFit_module.cc b/Mu2eKinKal/src/LoopHelixFit_module.cc index d44d92fcb8..e386deaa0a 100644 --- a/Mu2eKinKal/src/LoopHelixFit_module.cc +++ b/Mu2eKinKal/src/LoopHelixFit_module.cc @@ -108,7 +108,6 @@ namespace mu2e { using EXINGPTR = std::shared_ptr; using EXINGCOL = std::vector; - using KKMaterialConfig = KKMaterial::Config; using Name = fhicl::Name; using Comment = fhicl::Comment; @@ -125,7 +124,6 @@ namespace mu2e { fhicl::Table kkfitSettings { Name("KKFitSettings") }; fhicl::Table fitSettings { Name("FitSettings") }; fhicl::Table extSettings { Name("ExtensionSettings") }; - fhicl::Table matSettings { Name("MaterialSettings") }; fhicl::OptionalTable finalSettings { Name("FinalSettings") }; fhicl::OptionalTable extrapSettings { Name("ExtrapolationSettings") }; // LoopHelix module specific config @@ -167,7 +165,6 @@ namespace mu2e { TrkFitDirection fdir_; bool usePDGCharge_; // use the pdg particle charge: otherwise use the helicity and direction to determine the charge KKFIT kkfit_; // fit helper - KKMaterial kkmat_; // material helper DMAT seedcov_; // seed covariance matrix double mass_; // particle mass int PDGcharge_; // PDG particle charge @@ -198,7 +195,6 @@ namespace mu2e { useHelixSlope_(settings().slopeSigThreshold(slopeSigThreshold_)), usePDGCharge_(settings().pdgCharge()), kkfit_(settings().kkfitSettings()), - kkmat_(settings().matSettings()), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), exconfig_(Mu2eKinKal::makeConfig(settings().extSettings())), fixedfield_(false) @@ -224,7 +220,7 @@ namespace mu2e { kkbf_ = std::move(std::make_unique(VEC3(0.0,0.0,bz))); } // setup extrapolation - if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings(),kkmat_); + if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings()); // setup optional fit finalization; this just updates the internals, not the fit result itself if(settings().finalSettings()){ @@ -288,10 +284,9 @@ namespace mu2e { // check the input if(fdir.fitDirection() != TrkFitDirection::FitDirection::downstream && fdir.fitDirection() != TrkFitDirection::FitDirection::upstream) throw cet::exception("RECO") << "mu2e::LoopHelixFit: Unknown helix propagation direction " << fdir.name(); - - // Retrieve event information - // calo geom + // geom GeomHandle calo_h; + GeomHandle kkmat_h; // find current proditions auto const& strawresponse = strawResponse_h_.getPtr(event.id()); auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); @@ -336,7 +331,7 @@ namespace mu2e { strawhits.reserve(strawHitIdxs.size()); KKSTRAWXINGCOL strawxings; strawxings.reserve(strawHitIdxs.size()); - if(!kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_.strawMaterial(), pseedtraj, chcol, strawHitIdxs, strawhits, strawxings)) { + if(!kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_h->strawMaterial(), pseedtraj, chcol, strawHitIdxs, strawhits, strawxings)) { if(print_>0) printf("[LoopHelixFit::%s] Failed to create a track\n", __func__); return nullptr; } @@ -358,7 +353,7 @@ namespace mu2e { __func__, goodfit, ktrk->fitStatus().chisq_.probability(), ktrk->strawHits().size(), ktrk->caloHits().size()); // if we have an extension schedule, extend. if(goodfit && exconfig_.schedule().size() > 0) { - kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_.strawMaterial(), chcol, *calo_h, cc_H, *ktrk ); + kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H, *ktrk ); goodfit = goodFit(*ktrk,seedtraj); // if finalizing, apply that now. if(goodfit && fconfig_.schedule().size() > 0){ diff --git a/Mu2eKinKal/src/RegrowKinematicLine_module.cc b/Mu2eKinKal/src/RegrowKinematicLine_module.cc index 93d6e278d4..179331d13c 100644 --- a/Mu2eKinKal/src/RegrowKinematicLine_module.cc +++ b/Mu2eKinKal/src/RegrowKinematicLine_module.cc @@ -79,7 +79,6 @@ namespace mu2e { using Mu2eKinKal::KKFinalConfig; using KKFitConfig = Mu2eKinKal::KKFitConfig; using KKModuleConfig = Mu2eKinKal::KKModuleConfig; - using KKMaterialConfig = KKMaterial::Config; using SDIS = std::set; using Name = fhicl::Name; @@ -91,7 +90,6 @@ namespace mu2e { fhicl::Atom indexMap {Name("StrawDigiIndexMap"), Comment("Map between original and reduced ComboHits") }; fhicl::OptionalAtom kalSeedMCAssns {Name("KalSeedMCAssns"), Comment("Association to KalSeedMC. If set, regrown KalSeeds will be associated with the same KalSeedMC as the original") }; fhicl::Table kkfitSettings { Name("KKFitSettings") }; - fhicl::Table matSettings { Name("MaterialSettings") }; fhicl::Table fitSettings { Name("RefitSettings") }; fhicl::OptionalTable extrapSettings { Name("ExtrapolationSettings") }; }; @@ -130,8 +128,6 @@ namespace mu2e { using DOMAINPTR = std::shared_ptr; using DOMAINCOL = std::set; - using KKMaterialConfig = KKMaterial::Config; - explicit RegrowKinematicLine(const Parameters& settings); void beginRun(art::Run& run) override; void produce(art::Event& event) override; @@ -143,7 +139,6 @@ namespace mu2e { std::unique_ptr kkbf_; Config config_; // refit configuration object, containing the fit schedule KKFIT kkfit_; - KKMaterial kkmat_; art::ProductToken kseedcol_T_; art::ProductToken chcol_T_; art::ProductToken indexmap_T_; @@ -156,7 +151,6 @@ namespace mu2e { debug_(settings().debug()), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), kkfit_(settings().kkfitSettings()), - kkmat_(settings().matSettings()), kseedcol_T_(consumes(settings().kalSeedCollection())), chcol_T_(consumes(settings().comboHitCollection())), indexmap_T_(consumes(settings().indexMap())), @@ -164,7 +158,7 @@ namespace mu2e { { produces(); produces(); - if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings(),kkmat_); + if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings()); if( fillMCAssns_){ consumes(ksmca_T_); produces (); @@ -185,6 +179,8 @@ namespace mu2e { auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); GeomHandle nominalTracker_h; GeomHandle calo_h; + GeomHandle kkmat_h; + // find input event data auto kseed_H = event.getValidHandle(kseedcol_T_); const auto& kseedcol = *kseed_H; @@ -218,7 +214,7 @@ namespace mu2e { DOMAINCOL domains; // create the trajectory. This is done here to be strongly typed auto goodhits = kkfit_.regrowComponents(kseed, chcol, indexmap, - *tracker,*calo_h,*strawresponse,*kkbf_, kkmat_.strawMaterial(), + *tracker,*calo_h,*strawresponse,*kkbf_, kkmat_h->strawMaterial(), trajptr, strawhits, calohits, strawxings, domains); if(debug_ > 0){ std::cout << "Regrew " << strawhits.size() << " straw hits, " << strawxings.size() << " straw xings, " << calohits.size() << " CaloHits and " << domains.size() << " domains, status = " << goodhits << std::endl; diff --git a/Mu2eKinKal/src/RegrowLoopHelix_module.cc b/Mu2eKinKal/src/RegrowLoopHelix_module.cc index 2f648f2b6a..a4584dfe52 100644 --- a/Mu2eKinKal/src/RegrowLoopHelix_module.cc +++ b/Mu2eKinKal/src/RegrowLoopHelix_module.cc @@ -92,7 +92,6 @@ namespace mu2e { fhicl::Atom indexMap {Name("StrawDigiIndexMap"), Comment("Map between original and reduced ComboHits") }; fhicl::OptionalAtom kalSeedMCAssns {Name("KalSeedMCAssns"), Comment("Association to KalSeedMC. If set, regrown KalSeeds will be associated with the same KalSeedMC as the original") }; fhicl::Table kkfitSettings { Name("KKFitSettings") }; - fhicl::Table matSettings { Name("MaterialSettings") }; fhicl::Table fitSettings { Name("RefitSettings") }; fhicl::Atom extend {Name("Extend"), Comment("Extend the fit") }; @@ -146,7 +145,6 @@ namespace mu2e { std::unique_ptr kkbf_; Config config_; // refit configuration object, containing the fit schedule KKFIT kkfit_; - KKMaterial kkmat_; art::ProductToken kseedcol_T_; art::ProductToken chcol_T_; art::ProductToken cccol_T_; @@ -161,7 +159,6 @@ namespace mu2e { debug_(settings().debug()), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), kkfit_(settings().kkfitSettings()), - kkmat_(settings().matSettings()), kseedcol_T_(consumes(settings().kalSeedCollection())), chcol_T_(consumes(settings().comboHitCollection())), cccol_T_(mayConsume(settings().caloClusterCollection())), @@ -171,7 +168,7 @@ namespace mu2e { { produces(); produces(); - if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings(),kkmat_); + if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings()); if( fillMCAssns_){ consumes(ksmca_T_); produces (); @@ -192,6 +189,7 @@ namespace mu2e { auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); GeomHandle nominalTracker_h; GeomHandle calo_h; + GeomHandle kkmat_h; // find input event data auto kseed_H = event.getValidHandle(kseedcol_T_); const auto& kseedcol = *kseed_H; @@ -226,7 +224,7 @@ namespace mu2e { DOMAINCOL domains; // create the trajectory. This is done here to be strongly typed auto goodhits = kkfit_.regrowComponents(kseed, chcol, indexmap, - *tracker,*calo_h,*strawresponse,*kkbf_, kkmat_.strawMaterial(), + *tracker,*calo_h,*strawresponse,*kkbf_, kkmat_h->strawMaterial(), trajptr, strawhits, calohits, strawxings, domains); if(debug_ > 1){ std::cout << "Regrew " << strawhits.size() << " straw hits, " << strawxings.size() << " straw xings, " << calohits.size() << " CaloHits and " << domains.size() << " domains, status = " << goodhits << std::endl; @@ -245,7 +243,7 @@ namespace mu2e { auto ktrk = std::make_unique(config_,*kkbf_,kseed.particle(),trajptr,strawhits,strawxings,calohits,domains); if(ktrk && ktrk->fitStatus().usable()){ if(debug_ > 0) std::cout << "RegrowLoopHelix: successful track refit" << std::endl; - if(extend_)kkfit_.extendTrack(config_,*kkbf_, *tracker,*strawresponse, kkmat_.strawMaterial(), chcol, *calo_h, cc_H , *ktrk ); + if(extend_)kkfit_.extendTrack(config_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H , *ktrk ); if(ktrk->fitStatus().usable()){ // extrapolate as requested if(extrap_)extrap_->extrapolate(*ktrk); From 0004c25f1e9e5693fd09e8fe61f99f3298450bf9 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Sat, 16 May 2026 14:44:55 -0700 Subject: [PATCH 08/12] Fix --- Mu2eKinKal/src/KinematicLineFit_module.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Mu2eKinKal/src/KinematicLineFit_module.cc b/Mu2eKinKal/src/KinematicLineFit_module.cc index 60f135404c..a2143c706b 100644 --- a/Mu2eKinKal/src/KinematicLineFit_module.cc +++ b/Mu2eKinKal/src/KinematicLineFit_module.cc @@ -259,6 +259,8 @@ namespace mu2e { void KinematicLineFit::produce(art::Event& event ) { GeomHandle calo_h; GeomHandle nominalTracker_h; + GeomHandle kkmat_h; + // find current proditions auto const& strawresponse = strawResponse_h_.getPtr(event.id()); auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); From 9edc8a0bfdde16faa0d619677818cc4982dd1bc7 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Mon, 18 May 2026 11:14:39 -0700 Subject: [PATCH 09/12] Remove unneeded ProditionsEntity dependency. Cleanup fcl --- GeometryService/src/KinKalGeomMaker.cc | 16 ++++++++-------- KinKalGeom/fcl/prolog.fcl | 18 ++++++++++++++++++ KinKalGeom/inc/KKMaterial.hh | 4 +--- KinKalGeom/inc/KinKalGeom.hh | 5 ++--- KinKalGeom/src/KKMaterial.cc | 2 +- KinKalGeom/src/KinKalGeom.cc | 1 - Mu2eKinKal/fcl/prolog.fcl | 1 + TrackerConditions/fcl/prolog.fcl | 14 -------------- fcl/standardServices.fcl | 26 +++++++------------------- 9 files changed, 38 insertions(+), 49 deletions(-) create mode 100644 KinKalGeom/fcl/prolog.fcl diff --git a/GeometryService/src/KinKalGeomMaker.cc b/GeometryService/src/KinKalGeomMaker.cc index 9d33dbccdb..683b22a997 100644 --- a/GeometryService/src/KinKalGeomMaker.cc +++ b/GeometryService/src/KinKalGeomMaker.cc @@ -81,14 +81,14 @@ namespace mu2e { void KinKalGeomMaker::makeDS() { GeomHandle det; GeomHandle ds; - std::cout << "DS Cryo or " << ds->rOut1() << "," << ds->rOut2() << " ir " << ds->rIn1()<<","<< ds->rIn2() << " halfl " << ds->halfLength() - << " zpos " << ds->position().z() << " material " << ds->material() << std::endl; - std::cout << "DS shield or " << ds->shield_rOut1() << "," << ds->shield_rOut2() << " ir " << ds->shield_rIn1()<<","<< ds->shield_rIn2() << " halfl " << ds->shield_halfLength() << " material " << ds->shield_material() << std::endl; - std::cout << "DS ncoils " << ds->nCoils() << std::endl; - for(size_t icoil = 0; icoil < static_cast(ds->nCoils()); icoil++){ - std::cout << "DS coil ir " << ds->coil_rIn() << " or " << ds->coil_rOut()[icoil] << " length " << ds->coil_zLength()[icoil] << " zpos " << ds->coil_zPosition()[icoil] - << " material " << ds->coil_materials()[icoil] << std::endl; - } +// std::cout << "DS Cryo or " << ds->rOut1() << "," << ds->rOut2() << " ir " << ds->rIn1()<<","<< ds->rIn2() << " halfl " << ds->halfLength() +// << " zpos " << ds->position().z() << " material " << ds->material() << std::endl; +// std::cout << "DS shield or " << ds->shield_rOut1() << "," << ds->shield_rOut2() << " ir " << ds->shield_rIn1()<<","<< ds->shield_rIn2() << " halfl " << ds->shield_halfLength() << " material " << ds->shield_material() << std::endl; +// std::cout << "DS ncoils " << ds->nCoils() << std::endl; +// for(size_t icoil = 0; icoil < static_cast(ds->nCoils()); icoil++){ +// std::cout << "DS coil ir " << ds->coil_rIn() << " or " << ds->coil_rOut()[icoil] << " length " << ds->coil_zLength()[icoil] << " zpos " << ds->coil_zPosition()[icoil] +// << " material " << ds->coil_materials()[icoil] << std::endl; +// } //DS Cryo or 1303,1328 ir 950,969.05 halfl 5450 zpos 8689 material StainlessSteel //DS shield or 1237.3,1250 ir 1010,1022.7 halfl 5287.7 material G4_Al //DS ncoils 11 diff --git a/KinKalGeom/fcl/prolog.fcl b/KinKalGeom/fcl/prolog.fcl new file mode 100644 index 0000000000..ddf2db8258 --- /dev/null +++ b/KinKalGeom/fcl/prolog.fcl @@ -0,0 +1,18 @@ +BEGIN_PROLOG +KinKalGeom : { + KKMaterial : { + elements : "Offline/TrackerConditions/data/ElementsList.data" + isotopes : "Offline/TrackerConditions/data/IsotopesList.data" + materials : "Offline/TrackerConditions/data/MaterialsList.data" + strawGasMaterialName : "straw-gas" + strawWallMaterialName : "straw-wall" + strawWireMaterialName : "straw-wire" + IPAMaterialName : "HDPE" + STMaterialName : "Target" + IonizationEnergyLossMode : 1 # Moyal mean + SolidScatteringFraction : 0.999999 + GasScatteringFraction : 0.9999999 + ElectronBrehmsFraction : 0.04 + } +} +END_PROLOG diff --git a/KinKalGeom/inc/KKMaterial.hh b/KinKalGeom/inc/KKMaterial.hh index 13b83d62eb..ef8a221ea8 100644 --- a/KinKalGeom/inc/KKMaterial.hh +++ b/KinKalGeom/inc/KKMaterial.hh @@ -14,13 +14,12 @@ // mu2e #include "Offline/ConfigTools/inc/ConfigFileLookupPolicy.hh" #include "Offline/Mu2eInterfaces/inc/Detector.hh" -#include "Offline/Mu2eInterfaces/inc/ProditionsEntity.hh" #include #include namespace mu2e { - class KKMaterial : public MatEnv::FileFinderInterface, public Detector, public ProditionsEntity { + class KKMaterial : public MatEnv::FileFinderInterface, public Detector { public: using Name = fhicl::Name; using Comment = fhicl::Comment; @@ -29,7 +28,6 @@ namespace mu2e { fhicl::Atom isotopes { Name("isotopes"), Comment("Filename for istotopes information")}; fhicl::Atom elements { Name("elements"), Comment("Filename for elements information") }; fhicl::Atom materials { Name("materials"), Comment("Filename for materials information") }; - fhicl::Atom eloss { Name("ELossMode"), Comment("Energy Loss model (0=MPV, 1=Moyal"),MatEnv::DetMaterial::moyalmean }; fhicl::Atom strawGasMaterialName{ Name("strawGasMaterialName"), Comment("strawGasMaterialName") }; fhicl::Atom strawWallMaterialName{ Name("strawWallMaterialName"), Comment("strawWallMaterialName") }; fhicl::Atom strawWireMaterialName{ Name("strawWireMaterialName"), Comment("strawWireMaterialName") }; diff --git a/KinKalGeom/inc/KinKalGeom.hh b/KinKalGeom/inc/KinKalGeom.hh index 3d9a7e7dc2..87c04e790e 100644 --- a/KinKalGeom/inc/KinKalGeom.hh +++ b/KinKalGeom/inc/KinKalGeom.hh @@ -13,12 +13,11 @@ #include "Offline/DataProducts/inc/SurfaceId.hh" #include "KinKal/Geometry/Surface.hh" #include "Offline/Mu2eInterfaces/inc/Detector.hh" -#include "Offline/Mu2eInterfaces/inc/ProditionsEntity.hh" #include #include #include namespace mu2e { - class KinKalGeom : public Detector, public ProditionsEntity { + class KinKalGeom : public Detector { public: using SurfacePtr = std::shared_ptr; using SurfacePair =std::pair; @@ -26,7 +25,7 @@ namespace mu2e { using SurfacePairIter = std::multimap::const_iterator; using KKGMap = std::multimap; // default constructor, now using GeometryService to create content - KinKalGeom(); + KinKalGeom(){} // accessor to the raw map auto const& map() const { return map_; } // find all surfaces that match an Id. Return vector can have >1 entry if the wildcard index (-1) is provided diff --git a/KinKalGeom/src/KKMaterial.cc b/KinKalGeom/src/KKMaterial.cc index 8c390445d6..e80edefdb5 100644 --- a/KinKalGeom/src/KKMaterial.cc +++ b/KinKalGeom/src/KKMaterial.cc @@ -7,7 +7,7 @@ namespace mu2e { using MatDBInfo = MatEnv::MatDBInfo; using MatEnv::DetMaterial; - KKMaterial::KKMaterial(KKMaterial::Config const& matconfig) : ProditionsEntity("KKMaterial"), + KKMaterial::KKMaterial(KKMaterial::Config const& matconfig) : elementsBaseName_(matconfig.elements()), isotopesBaseName_(matconfig.isotopes()), materialsBaseName_(matconfig.materials()), diff --git a/KinKalGeom/src/KinKalGeom.cc b/KinKalGeom/src/KinKalGeom.cc index e6c21fd2cd..13bf8adcff 100644 --- a/KinKalGeom/src/KinKalGeom.cc +++ b/KinKalGeom/src/KinKalGeom.cc @@ -1,7 +1,6 @@ #include "Offline/KinKalGeom/inc/KinKalGeom.hh" #include "cetlib_except/exception.h" namespace mu2e { - KinKalGeom::KinKalGeom() : ProditionsEntity("KinKalGeom") {} void KinKalGeom::surfaces(SurfaceId const& id, SurfacePairCollection& surfs) const { // find all surfaces that match the input, including index wildcard diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 808ff07c56..03fc9879a2 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -1,6 +1,7 @@ #include "Offline/TrkReco/fcl/Particle.fcl" BEGIN_PROLOG Mu2eKinKal : { + # # general configuration # KKFIT: { diff --git a/TrackerConditions/fcl/prolog.fcl b/TrackerConditions/fcl/prolog.fcl index edc76089bd..2fb1239208 100644 --- a/TrackerConditions/fcl/prolog.fcl +++ b/TrackerConditions/fcl/prolog.fcl @@ -1,18 +1,4 @@ BEGIN_PROLOG -Mu2eMaterial : { - verbose : 0 - # Location of dictionary files for the material model. - elements : "Offline/TrackerConditions/data/ElementsList.data" - isotopes : "Offline/TrackerConditions/data/IsotopesList.data" - materials : "Offline/TrackerConditions/data/MaterialsList.data" - strawGasMaterialName : "straw-gas" - strawWallMaterialName : "straw-wall" - strawWireMaterialName : "straw-wire" - dahlLynchScatteringFraction : 0.995 - intersectionTolerance : 0.001 - strawElementOffset : 0.25 - maximumIntersectionRadiusFraction : 0.96 -} Mu2eDetector : { verbose : 0 diff --git a/fcl/standardServices.fcl b/fcl/standardServices.fcl index 5f121f475e..bea908f6c4 100644 --- a/fcl/standardServices.fcl +++ b/fcl/standardServices.fcl @@ -5,6 +5,7 @@ #include "Offline/fcl/minimalMessageService.fcl" #include "Offline/DbService/fcl/prolog.fcl" #include "Offline/ProditionsService/fcl/prolog.fcl" +#include "Offline/KinKalGeom/fcl/prolog.fcl" BEGIN_PROLOG @@ -37,25 +38,12 @@ Services : { # define services for specific tasks: these are components needed for # a complete job Core : { - message : @local::default_message - GeometryService : { - inputFile : "Offline/Mu2eG4/geom/geom_common.txt" - bFieldFile : "Offline/Mu2eG4/geom/bfgeom_v01.txt" - simulatedDetector : { tool_type: "Mu2e" } - KinKalMaterial : { - elements : "Offline/TrackerConditions/data/ElementsList.data" - isotopes : "Offline/TrackerConditions/data/IsotopesList.data" - materials : "Offline/TrackerConditions/data/MaterialsList.data" - strawGasMaterialName : "straw-gas" - strawWallMaterialName : "straw-wall" - strawWireMaterialName : "straw-wire" - IPAMaterialName : "HDPE" - STMaterialName : "Target" - IonizationEnergyLossMode : 1 # Moyal mean - SolidScatteringFraction : 0.999999 - GasScatteringFraction : 0.9999999 - ElectronBrehmsFraction : 0.04 - } + message : @local::default_message + GeometryService : { + inputFile : "Offline/Mu2eG4/geom/geom_common.txt" + bFieldFile : "Offline/Mu2eG4/geom/bfgeom_v01.txt" + simulatedDetector : { tool_type: "Mu2e" } + KinKalMaterial : @local::KinKalGeom.KKMaterial } GlobalConstantsService : { inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt" } DbService : @local::DbEmpty From a4d9616c06bd16ff35d8c87cc290cd1fb178adda Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Mon, 18 May 2026 12:06:49 -0700 Subject: [PATCH 10/12] Cleanup interface --- Mu2eKinKal/inc/KKFit.hh | 35 ++++++++++++-------- Mu2eKinKal/src/CentralHelixFit_module.cc | 6 ++-- Mu2eKinKal/src/KinematicLineFit_module.cc | 5 ++- Mu2eKinKal/src/LoopHelixFit_module.cc | 6 ++-- Mu2eKinKal/src/RegrowKinematicLine_module.cc | 4 +-- Mu2eKinKal/src/RegrowLoopHelix_module.cc | 9 ++--- 6 files changed, 31 insertions(+), 34 deletions(-) diff --git a/Mu2eKinKal/inc/KKFit.hh b/Mu2eKinKal/inc/KKFit.hh index 32e47d6aba..36f5d2bb87 100644 --- a/Mu2eKinKal/inc/KKFit.hh +++ b/Mu2eKinKal/inc/KKFit.hh @@ -7,7 +7,7 @@ #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" -#include "Offline/KinKalGeom/inc/KKStrawMaterial.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKCaloHit.hh" #include "Offline/Mu2eKinKal/inc/KKFitUtilities.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" @@ -96,18 +96,18 @@ namespace mu2e { explicit KKFit(KKFitConfig const& fitconfig); // helper functions used to create components of the fit // Make KKStrawHits from ComboHits - bool makeStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf, KKStrawMaterial const& smat, + bool makeStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf, PTRAJ const& ptraj, ComboHitCollection const& chcol, StrawHitIndexCollection const& strawHitIdxs, KKSTRAWHITCOL& hits, KKSTRAWXINGCOL& exings) const; // regrow KKTrack components from a KalSeed bool regrowComponents(KalSeed const& kseed, ComboHitCollection const& chcol, mu2e::IndexMap const& strawindexmap, - Tracker const& tracker,Calorimeter const& calo, StrawResponse const& strawresponse,BFieldMap const& kkbf, KKStrawMaterial const& smat, + Tracker const& tracker,Calorimeter const& calo, StrawResponse const& strawresponse,BFieldMap const& kkbf, PTRAJPTR& ptraj, KKSTRAWHITCOL& strawhits, KKCALOHITCOL& calohits, KKSTRAWXINGCOL& exings, DOMAINCOL& domains) const; std::shared_ptr caloAxis(CaloCluster const& cluster, Calorimeter const& calo) const; // should come from CaloCluster TODO bool makeCaloHit(CCPtr const& cluster, Calorimeter const& calo, PTRAJ const& pktraj, KKCALOHITCOL& hits) const; // extend a track with a new configuration, optionally searching for and adding hits and straw material void extendTrack(Config const& config, BFieldMap const& kkbf, Tracker const& tracker, - StrawResponse const& strawresponse, KKStrawMaterial const& smat, ComboHitCollection const& chcol, + StrawResponse const& strawresponse, ComboHitCollection const& chcol, Calorimeter const& calo, CCHandle const& cchandle, KKTRK& kktrk) const; // sample the fit at the specificed surfaces @@ -122,9 +122,9 @@ namespace mu2e { auto const& strawHitClusterer() const { return shclusterer_; } private: void fillTrackerInfo(Tracker const& tracker) const; - void addStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf, KKStrawMaterial const& smat, + void addStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf, KKTRK const& kktrk, ComboHitCollection const& chcol, KKSTRAWHITCOL& hits,KKSTRAWXINGCOL& addexings) const; - void addStraws(Tracker const& tracker, KKStrawMaterial const& smat, KKTRK const& kktrk, KKSTRAWXINGCOL& addexings) const; + void addStraws(Tracker const& tracker, KKTRK const& kktrk, KKSTRAWXINGCOL& addexings) const; void addCaloHit(Calorimeter const& calo, KKTRK& kktrk, CCHandle cchandle, KKCALOHITCOL& hits) const; void sampleFit(KKTRK const& kktrk,KalIntersectionCollection& inters) const; // sample fit at the surfaces specified in the config int print_; @@ -212,10 +212,12 @@ namespace mu2e { } } - template bool KKFit::makeStrawHits(Tracker const& tracker,StrawResponse const& strawresponse,BFieldMap const& kkbf, KKStrawMaterial const& smat, + template bool KKFit::makeStrawHits(Tracker const& tracker,StrawResponse const& strawresponse,BFieldMap const& kkbf, PTRAJ const& ptraj, ComboHitCollection const& chcol, StrawHitIndexCollection const& strawHitIdxs, KKSTRAWHITCOL& hits, KKSTRAWXINGCOL& exings) const { unsigned ngood(0); + GeomHandle kkmat_h; + auto const& smat = kkmat_h->strawMaterial(); // loop over the individual straw combo hits for(auto strawidx : strawHitIdxs) { const ComboHit& combohit(chcol.at(strawidx)); @@ -252,8 +254,10 @@ namespace mu2e { template bool KKFit::regrowComponents(KalSeed const& kseed, // primary event input ComboHitCollection const& chcol, mu2e::IndexMap const& strawindexmap, // ancillary event input Tracker const& tracker,Calorimeter const& calo, // geometries - StrawResponse const& strawresponse,BFieldMap const& kkbf, KKStrawMaterial const& smat, // other conditions + StrawResponse const& strawresponse,BFieldMap const& kkbf, // other conditions PTRAJPTR& ptraj, KKSTRAWHITCOL& strawhits, KKCALOHITCOL& calohits, KKSTRAWXINGCOL& exings, DOMAINCOL& domains) const { // return values + GeomHandle kkmat_h; + auto const& smat = kkmat_h->strawMaterial(); unsigned ngood(0), nactive(0), nsactive(0); // loop over the TrkStrawHitSeeds in this KalSeed for(auto const& tshs : kseed.hits()) { @@ -349,14 +353,14 @@ namespace mu2e { } template void KKFit::extendTrack(Config const& exconfig, BFieldMap const& kkbf, Tracker const& tracker, - StrawResponse const& strawresponse, KKStrawMaterial const& smat, ComboHitCollection const& chcol, + StrawResponse const& strawresponse, ComboHitCollection const& chcol, Calorimeter const& calo, CCHandle const& cchandle, KKTRK& kktrk) const { KKSTRAWHITCOL addstrawhits; KKCALOHITCOL addcalohits; KKSTRAWXINGCOL addstrawxings; - if(addhits_)addStrawHits(tracker, strawresponse, kkbf, smat, kktrk, chcol, addstrawhits, addstrawxings ); - if(matcorr_ && addmat_)addStraws(tracker, smat, kktrk, addstrawxings); + if(addhits_)addStrawHits(tracker, strawresponse, kkbf, kktrk, chcol, addstrawhits, addstrawxings ); + if(matcorr_ && addmat_)addStraws(tracker, kktrk, addstrawxings); if(addhits_ && usecalo_ && kktrk.caloHits().size()==0)addCaloHit(calo, kktrk, cchandle, addcalohits); if(print_ > 1){ std::cout << "KKTrk extension adding " @@ -367,8 +371,10 @@ namespace mu2e { kktrk.extendTrack(exconfig,addstrawhits,addstrawxings,addcalohits); } - template void KKFit::addStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf, KKStrawMaterial const& smat, + template void KKFit::addStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf, KKTRK const& kktrk, ComboHitCollection const& chcol, KKSTRAWHITCOL& addhits, KKSTRAWXINGCOL& addexings) const { + GeomHandle kkmat_h; + auto const& smat = kkmat_h->strawMaterial(); auto const& ptraj = kktrk.fitTraj(); // add the buffer to the time range; this defines the search range for new hits TimeRange brange(ptraj.range().begin()-maxStrawHitDt_, ptraj.range().end()+maxStrawHitDt_); @@ -429,8 +435,11 @@ namespace mu2e { } } - template void KKFit::addStraws(Tracker const& tracker, KKStrawMaterial const& smat, KKTRK const& kktrk, + template void KKFit::addStraws(Tracker const& tracker, KKTRK const& kktrk, KKSTRAWXINGCOL& addexings) const { + GeomHandle kkmat_h; + auto const& smat = kkmat_h->strawMaterial(); + // this algorithm assumes the track never hits the same straw twice. That could be violated by reflecting tracks, and could be addressed // by including the time of the Xing as part of its identity. That would slow things down so it remains to be proven it's a problem TODO // build the set of existing straws diff --git a/Mu2eKinKal/src/CentralHelixFit_module.cc b/Mu2eKinKal/src/CentralHelixFit_module.cc index f6cef066ed..90ba0e1ef2 100644 --- a/Mu2eKinKal/src/CentralHelixFit_module.cc +++ b/Mu2eKinKal/src/CentralHelixFit_module.cc @@ -52,7 +52,6 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" -#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" @@ -242,7 +241,6 @@ namespace mu2e { void CentralHelixFit::produce(art::Event& event ) { GeomHandle calo_h; GeomHandle nominalTracker_h; - GeomHandle kkmat_h; // find current proditions auto const& strawresponse = strawResponse_h_.getPtr(event.id()); auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); @@ -316,7 +314,7 @@ namespace mu2e { strawhits.reserve(strawHitIdxs.size()); KKSTRAWXINGCOL strawxings; strawxings.reserve(strawHitIdxs.size()); - kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_h->strawMaterial(), pseedtraj, chcol, strawHitIdxs, strawhits, strawxings); + kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, pseedtraj, chcol, strawHitIdxs, strawhits, strawxings); // optionally (and if present) add the CaloCluster as a constraint // verify the cluster looks physically reasonable before adding it TODO! Or, let the KKCaloHit updater do it TODO KKCALOHITCOL calohits; @@ -332,7 +330,7 @@ namespace mu2e { // if we have an extension schedule, extend. if(goodfit && exconfig_.schedule().size() > 0) { // std::cout << "EXTENDING TRACK " << event.id() << " " << index << std::endl; - kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H, *kktrk ); + kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, chcol, *calo_h, cc_H, *kktrk ); goodfit = goodFit(*kktrk); } diff --git a/Mu2eKinKal/src/KinematicLineFit_module.cc b/Mu2eKinKal/src/KinematicLineFit_module.cc index 9f309c0ba8..28c96e40e3 100644 --- a/Mu2eKinKal/src/KinematicLineFit_module.cc +++ b/Mu2eKinKal/src/KinematicLineFit_module.cc @@ -259,7 +259,6 @@ namespace mu2e { void KinematicLineFit::produce(art::Event& event ) { GeomHandle calo_h; GeomHandle nominalTracker_h; - GeomHandle kkmat_h; // find current proditions auto const& strawresponse = strawResponse_h_.getPtr(event.id()); @@ -300,7 +299,7 @@ namespace mu2e { KKSTRAWXINGCOL strawxings; strawhits.reserve(strawHitIdxs.size()); strawxings.reserve(strawHitIdxs.size()); - kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_h->strawMaterial(), pseedtraj, *chcolptr, strawHitIdxs, strawhits, strawxings); + kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, pseedtraj, *chcolptr, strawHitIdxs, strawhits, strawxings); //here KKCALOHITCOL calohits; @@ -321,7 +320,7 @@ namespace mu2e { auto kktrk = make_unique(config_,*kkbf_,seedtraj,fpart_,kkfit_.strawHitClusterer(),strawhits,strawxings,calohits,paramconstraints_); auto goodfit = goodFit(*kktrk); if(goodfit && exconfig_.schedule().size() > 0){ - kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H, *kktrk ); + kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, chcol, *calo_h, cc_H, *kktrk ); } goodfit = goodFit(*kktrk); // extrapolate as required diff --git a/Mu2eKinKal/src/LoopHelixFit_module.cc b/Mu2eKinKal/src/LoopHelixFit_module.cc index fca260082a..1ecec5900f 100644 --- a/Mu2eKinKal/src/LoopHelixFit_module.cc +++ b/Mu2eKinKal/src/LoopHelixFit_module.cc @@ -52,7 +52,6 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" -#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" @@ -286,7 +285,6 @@ namespace mu2e { throw cet::exception("RECO") << "mu2e::LoopHelixFit: Unknown helix propagation direction " << fdir.name(); // geom GeomHandle calo_h; - GeomHandle kkmat_h; // find current proditions auto const& strawresponse = strawResponse_h_.getPtr(event.id()); auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); @@ -331,7 +329,7 @@ namespace mu2e { strawhits.reserve(strawHitIdxs.size()); KKSTRAWXINGCOL strawxings; strawxings.reserve(strawHitIdxs.size()); - if(!kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_h->strawMaterial(), pseedtraj, chcol, strawHitIdxs, strawhits, strawxings)) { + if(!kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, pseedtraj, chcol, strawHitIdxs, strawhits, strawxings)) { if(print_>0) printf("[LoopHelixFit::%s] Failed to create a track\n", __func__); return nullptr; } @@ -353,7 +351,7 @@ namespace mu2e { __func__, goodfit, ktrk->fitStatus().chisq_.probability(), ktrk->strawHits().size(), ktrk->caloHits().size()); // if we have an extension schedule, extend. if(goodfit && exconfig_.schedule().size() > 0) { - kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H, *ktrk ); + kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, chcol, *calo_h, cc_H, *ktrk ); goodfit = goodFit(*ktrk,seedtraj); // if finalizing, apply that now. if(goodfit && fconfig_.schedule().size() > 0){ diff --git a/Mu2eKinKal/src/RegrowKinematicLine_module.cc b/Mu2eKinKal/src/RegrowKinematicLine_module.cc index 179331d13c..ff02d69854 100644 --- a/Mu2eKinKal/src/RegrowKinematicLine_module.cc +++ b/Mu2eKinKal/src/RegrowKinematicLine_module.cc @@ -54,7 +54,6 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" -#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" @@ -179,7 +178,6 @@ namespace mu2e { auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); GeomHandle nominalTracker_h; GeomHandle calo_h; - GeomHandle kkmat_h; // find input event data auto kseed_H = event.getValidHandle(kseedcol_T_); @@ -214,7 +212,7 @@ namespace mu2e { DOMAINCOL domains; // create the trajectory. This is done here to be strongly typed auto goodhits = kkfit_.regrowComponents(kseed, chcol, indexmap, - *tracker,*calo_h,*strawresponse,*kkbf_, kkmat_h->strawMaterial(), + *tracker,*calo_h,*strawresponse,*kkbf_, trajptr, strawhits, calohits, strawxings, domains); if(debug_ > 0){ std::cout << "Regrew " << strawhits.size() << " straw hits, " << strawxings.size() << " straw xings, " << calohits.size() << " CaloHits and " << domains.size() << " domains, status = " << goodhits << std::endl; diff --git a/Mu2eKinKal/src/RegrowLoopHelix_module.cc b/Mu2eKinKal/src/RegrowLoopHelix_module.cc index a4584dfe52..085ed2166d 100644 --- a/Mu2eKinKal/src/RegrowLoopHelix_module.cc +++ b/Mu2eKinKal/src/RegrowLoopHelix_module.cc @@ -54,7 +54,6 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" -#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" @@ -79,7 +78,6 @@ namespace mu2e { using Mu2eKinKal::KKFinalConfig; using KKFitConfig = Mu2eKinKal::KKFitConfig; using KKModuleConfig = Mu2eKinKal::KKModuleConfig; - using KKMaterialConfig = KKMaterial::Config; using SDIS = std::set; using Name = fhicl::Name; @@ -132,8 +130,6 @@ namespace mu2e { using DOMAINPTR = std::shared_ptr; using DOMAINCOL = std::set; - using KKMaterialConfig = KKMaterial::Config; - explicit RegrowLoopHelix(const Parameters& settings); void beginRun(art::Run& run) override; void produce(art::Event& event) override; @@ -189,7 +185,6 @@ namespace mu2e { auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); GeomHandle nominalTracker_h; GeomHandle calo_h; - GeomHandle kkmat_h; // find input event data auto kseed_H = event.getValidHandle(kseedcol_T_); const auto& kseedcol = *kseed_H; @@ -224,7 +219,7 @@ namespace mu2e { DOMAINCOL domains; // create the trajectory. This is done here to be strongly typed auto goodhits = kkfit_.regrowComponents(kseed, chcol, indexmap, - *tracker,*calo_h,*strawresponse,*kkbf_, kkmat_h->strawMaterial(), + *tracker,*calo_h,*strawresponse,*kkbf_, trajptr, strawhits, calohits, strawxings, domains); if(debug_ > 1){ std::cout << "Regrew " << strawhits.size() << " straw hits, " << strawxings.size() << " straw xings, " << calohits.size() << " CaloHits and " << domains.size() << " domains, status = " << goodhits << std::endl; @@ -243,7 +238,7 @@ namespace mu2e { auto ktrk = std::make_unique(config_,*kkbf_,kseed.particle(),trajptr,strawhits,strawxings,calohits,domains); if(ktrk && ktrk->fitStatus().usable()){ if(debug_ > 0) std::cout << "RegrowLoopHelix: successful track refit" << std::endl; - if(extend_)kkfit_.extendTrack(config_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H , *ktrk ); + if(extend_)kkfit_.extendTrack(config_,*kkbf_, *tracker,*strawresponse, chcol, *calo_h, cc_H , *ktrk ); if(ktrk->fitStatus().usable()){ // extrapolate as requested if(extrap_)extrap_->extrapolate(*ktrk); From aef8ba89f7ac04c529cb34efd4db542db1946b8a Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Mon, 18 May 2026 13:50:28 -0700 Subject: [PATCH 11/12] Move toCRV to KKExtrap --- Mu2eKinKal/inc/KKExtrap.hh | 46 ++++++++++++++- Mu2eKinKal/inc/KKFitSettings.hh | 2 + Mu2eKinKal/src/KinematicLineFit_module.cc | 69 ++--------------------- 3 files changed, 52 insertions(+), 65 deletions(-) diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index f46290275d..cc430349fc 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -38,14 +38,17 @@ namespace mu2e { template bool extrapolateST(KKTrack& ktrk,TimeDir trkdir) const; template bool extrapolateTracker(KKTrack& ktrk,TimeDir tdir) const; template bool extrapolateTSDA(KKTrack& ktrk,TimeDir tdir) const; + template bool extrapolateTCRV(KKTrack& ktrk,TimeDir tdir) const; template void toOPA(KKTrack& ktrk, double tstart, TimeDir tdir) const; + template void toCRV(KKTrack& ktrk, TimeDir tdir) const; private: int debug_; - double btol_, intertol_, maxdt_; - bool backToTracker_, toOPA_, toTrackerEnds_, upstream_; + double btol_, intertol_, maxdt_, minv_; + bool backToTracker_, toOPA_, toTrackerEnds_, upstream_, toCRV_; double ipathick_ = 0.511; // ipa thickness: should come from geometry service TODO double stthick_ = 0.1056; // st foil thickness: should come from geometry service TODO + double tcrvthick_ = 150.0; // test CRV sector thickness: should come from geometry service TODO }; KKExtrap::KKExtrap(KKExtrapConfig const& extrapconfig) : @@ -53,8 +56,10 @@ namespace mu2e { btol_(extrapconfig.btol()), intertol_(extrapconfig.interTol()), maxdt_(extrapconfig.MaxDt()), + minv_(extrapconfig.MinV()), backToTracker_(extrapconfig.BackToTracker()), toOPA_(extrapconfig.ToOPA()), + toCRV_(extrapconfig.ToCRV()), toTrackerEnds_(extrapconfig.ToTrackerEnds()), upstream_(extrapconfig.Upstream()) {} @@ -95,6 +100,8 @@ namespace mu2e { } // optionally test for intersection with the OPA if(toOPA_)toOPA(ktrk,starttime,tdir); + // optionally test for intersection with the CRV + if(toCRV_)toCRV(ktrk,starttime,tdir); } } @@ -271,5 +278,40 @@ namespace mu2e { ktrk.addIntersection(OPASID,opainter); } } + + template void KKExtrap::toTCRV(KKTrack& ktrk,TimeDir tdir) const { + GeomHandle kkg_h; + GeomHandle kkmat_h; + // extrapolate to the extracted CRV + auto TCRV = ExtrapolateTCRV(maxdt_,btol_,intertol_,minv_,*kkg_h->TCRV(),debug_); + + auto const& ftraj = ktrk.fitTraj(); + static const SurfaceId TCRVSID("TCRV"); + auto dir0 = ftraj.direction(ftraj.t0()); + TimeDir tdir = (dir0.Y() > 0) ? TimeDir::forwards : TimeDir::backwards; + double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); + bool hadintersection = false; + do { + // iterate until the extrapolation condition is met + double time = starttime; + double tstart = time; + // cosmics are always (?) downwards, we can simplify this logic using that TODO + while(fabs(time-tstart) < TCRV.maxDt() && TCRV.needsExtrapolation(ftraj,tdir) ){ + TimeRange range = tdir == TimeDir::forwards ? TimeRange(time,time+TCRV.step()) : TimeRange(time-TCRV.step(),time); + ktrk.extendTraj(range); + time = tdir == TimeDir::forwards ? range.end() : range.begin(); + } + hadintersection = false; + if (TCRV.intersection().good()){ + hadintersection = true; + // we have a good intersection. Use this to create a Shell material Xing + auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); + // TODO add DS and shielding material + KKCRVXINGPTR crvxingptr = std::make_shared(TCRV.module(), TCRVSID, *kkmat_h->STMaterial(),TCRV.intersection(),reftrajptr,tcrvthick_,TCRV.interTolerance()); + ktrk.addTCRVXing(crvxingptr,tdir); + } + } while(hadintersection); + } + } #endif diff --git a/Mu2eKinKal/inc/KKFitSettings.hh b/Mu2eKinKal/inc/KKFitSettings.hh index 3a881c64f6..d1dd609765 100644 --- a/Mu2eKinKal/inc/KKFitSettings.hh +++ b/Mu2eKinKal/inc/KKFitSettings.hh @@ -112,10 +112,12 @@ namespace mu2e { fhicl::Atom btol { Name("BCorrTolerance"), Comment("Tolerance on BField correction momentum fractional accuracy (dimensionless)") }; fhicl::Atom interTol { Name("IntersectionTolerance"), Comment("Tolerance for surface intersections (mm)") }; fhicl::Atom MaxDt { Name("MaxDt"), Comment("Maximum time to extrapolate a fit (ns)") }; + fhicl::Atom MinV { Name("MinV"), Comment("Minimum Y vel to extrapolate a fit") }; fhicl::Atom BackToTracker { Name("BackToTracker"), Comment("Extrapolate reflecting tracks back to the tracker") }; fhicl::Atom ToTrackerEnds { Name("ToTrackerEnds"), Comment("Extrapolate tracks to the tracker ends") }; fhicl::Atom Upstream { Name("Upstream"), Comment("Extrapolate tracks upstream") }; fhicl::Atom ToOPA { Name("ToOPA"), Comment("Test tracks for intersection with the OPA") }; + fhicl::Atom ToCRV { Name("ToCRV"), Comment("Extrapolate tracks to the CRV") }; }; } } diff --git a/Mu2eKinKal/src/KinematicLineFit_module.cc b/Mu2eKinKal/src/KinematicLineFit_module.cc index 28c96e40e3..acf1594f6c 100644 --- a/Mu2eKinKal/src/KinematicLineFit_module.cc +++ b/Mu2eKinKal/src/KinematicLineFit_module.cc @@ -52,7 +52,7 @@ #include "Offline/Mu2eKinKal/inc/KKBField.hh" #include "Offline/Mu2eKinKal/inc/KKFitUtilities.hh" #include "Offline/Mu2eKinKal/inc/KKShellXing.hh" -#include "Offline/Mu2eKinKal/inc/ExtrapolateTCRV.hh" +#include "Offline/Mu2eKinKal/inc/KKExtrap.hh" // root #include "TH1F.h" #include "TTree.h" @@ -119,20 +119,12 @@ namespace mu2e { fhicl::Atom sampleTBuff { Name("SampleTimeBuffer"), Comment("Time buffer for sample intersections (nsec)") }; }; - // Extrapolation configuration - struct KKExtrapConfig { - fhicl::Atom MaxDt { Name("MaxDt"), Comment("Maximum time to extrapolate a fit") }; - fhicl::Atom MinV { Name("MinV"), Comment("Minimum Y vel to extrapolate a fit") }; - fhicl::Atom ToCRV { Name("ToCRV"), Comment("Extrapolate tracks to the CRV") }; - fhicl::Atom Debug { Name("Debug"), Comment("Debug level"), 0 }; - }; - struct GlobalConfig { fhicl::Table modSettings { Name("ModuleSettings") }; fhicl::Table mu2eSettings { Name("KKFitSettings") }; fhicl::Table fitSettings { Name("FitSettings") }; fhicl::Table extSettings { Name("ExtensionSettings") }; - fhicl::OptionalTable Extrapolation { Name("Extrapolation") }; + fhicl::OptionalTable Extrapolation { Name("ExtrapolationSettings") }; fhicl::OptionalAtom fitDirection { Name("FitDirection"), Comment("Particle direction to fit, either \"upstream\" or \"downstream\"")}; }; @@ -169,12 +161,9 @@ namespace mu2e { double intertol_; // surface intersection tolerance (mm) double sampletbuff_; // simple time buffer; replace this with extrapolation TODO bool sampleinrange_, sampleinbounds_; // require samples to be in range or on surface - bool extrapolate_, toCRV_; - double maxdt_ = 0.0, btol_ = 0.0, minv_ = 0.0; + std::unique_ptr extrap_; // extrapolation helper SurfaceIdCollection ssids_; KinKalGeom::SurfacePairCollection surfacess_to_sample_; // surfaces to sample the fit - int extrapdebug_ = 0; - double tcrvthick_ = 150.0; // CRV sector thickness: should come from geometry service TODO Config config_; // initial fit configuration object Config exconfig_; // extension configuration object }; @@ -192,7 +181,6 @@ namespace mu2e { sampletbuff_(settings().modSettings().sampleTBuff()), sampleinrange_(settings().modSettings().sampleInRange()), sampleinbounds_(settings().modSettings().sampleInBounds()), - extrapolate_(false), toCRV_(false), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), exconfig_(Mu2eKinKal::makeConfig(settings().extSettings())) { @@ -223,21 +211,9 @@ namespace mu2e { for(auto const& sidname : settings().modSettings().sampleSurfaces()) { ssids_.push_back(SurfaceId(sidname,-1)); // match all elements } - // configure extrapolation - if(settings().Extrapolation()){ - extrapolate_ = true; - toCRV_ = settings().Extrapolation()->ToCRV(); - // global configs - maxdt_ = settings().Extrapolation()->MaxDt(); - btol_ = settings().extSettings().btol(); // use the same BField cor. tolerance as in fit extension - minv_ = settings().Extrapolation()->MinV(); - extrapdebug_ = settings().Extrapolation()->Debug(); - } - - + // setup extrapolation + if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings()); if(print_ > 0) std::cout << config_; - - } KinematicLineFit::~KinematicLineFit(){} @@ -324,7 +300,7 @@ namespace mu2e { } goodfit = goodFit(*kktrk); // extrapolate as required - if(goodfit && extrapolate_) extrapolate(*kktrk); + if(goodfit && extrap_) extrap_->extrapolate(*kktrk); bool save = goodFit(*kktrk); if(save || saveall_){ TrkFitFlag fitflag(hptr->status()); @@ -404,38 +380,5 @@ namespace mu2e { } } } - - void KinematicLineFit::extrapolate(KKTRK& ktrk) const { - GeomHandle kkg_h; - GeomHandle kkmat_h; - // extrapolate to the extracted CRV. This function should be migrated to KKExtrap TODO - auto TCRV = ExtrapolateTCRV(maxdt_,btol_,intertol_,minv_,*kkg_h->TCRV(),extrapdebug_); - - auto const& ftraj = ktrk.fitTraj(); - static const SurfaceId TCRVSID("TCRV"); - auto dir0 = ftraj.direction(ftraj.t0()); - TimeDir tdir = (dir0.Y() > 0) ? TimeDir::forwards : TimeDir::backwards; - double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); - bool hadintersection = false; - do { - // iterate until the extrapolation condition is met - double time = starttime; - double tstart = time; - while(fabs(time-tstart) < TCRV.maxDt() && TCRV.needsExtrapolation(ftraj,tdir) ){ - TimeRange range = tdir == TimeDir::forwards ? TimeRange(time,time+TCRV.step()) : TimeRange(time-TCRV.step(),time); - ktrk.extendTraj(range); - time = tdir == TimeDir::forwards ? range.end() : range.begin(); - } - hadintersection = false; - if (TCRV.intersection().good()){ - hadintersection = true; - // we have a good intersection. Use this to create a Shell material Xing - auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); - // TODO add DS and shielding material - KKCRVXINGPTR crvxingptr = std::make_shared(TCRV.module(), TCRVSID, *kkmat_h->STMaterial(),TCRV.intersection(),reftrajptr,tcrvthick_,TCRV.interTolerance()); - ktrk.addTCRVXing(crvxingptr,tdir); - } - } while(hadintersection); - } } DEFINE_ART_MODULE(mu2e::KinematicLineFit) From b0d8ab61c06a90e30c07d780b12df6853aef7531 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Thu, 21 May 2026 07:52:38 -0700 Subject: [PATCH 12/12] resolve merge conflicts --- CalPatRec/inc/CalHelixFinderData.hh | 24 +- CalPatRec/inc/CalHelixFinder_module.hh | 6 +- CalPatRec/inc/CalHelixFinder_types.hh | 1 - CalPatRec/src/CalHelixFinderAlg.cc | 106 ++- CalPatRec/src/CalHelixFinderData.cc | 6 + CalPatRec/src/CalHelixFinder_module.cc | 685 +++++++++---------- Mu2eKinKal/inc/KKExtrap.hh | 46 +- Mu2eKinKal/inc/KKFit.hh | 35 +- Mu2eKinKal/inc/KKFitSettings.hh | 2 - Mu2eKinKal/src/CentralHelixFit_module.cc | 6 +- Mu2eKinKal/src/KinematicLineFit_module.cc | 74 +- Mu2eKinKal/src/LoopHelixFit_module.cc | 6 +- Mu2eKinKal/src/RegrowKinematicLine_module.cc | 4 +- Mu2eKinKal/src/RegrowLoopHelix_module.cc | 9 +- 14 files changed, 480 insertions(+), 530 deletions(-) diff --git a/CalPatRec/inc/CalHelixFinderData.hh b/CalPatRec/inc/CalHelixFinderData.hh index 858093d042..874bffc5bd 100644 --- a/CalPatRec/inc/CalHelixFinderData.hh +++ b/CalPatRec/inc/CalHelixFinderData.hh @@ -132,33 +132,27 @@ namespace mu2e { //----------------------------------------------------------------------------- // circle parameters; the z center is ignored. //----------------------------------------------------------------------------- - ::LsqSums4 _sxy; - ::LsqSums4 _szphi; + ::LsqSums4 _sxy; // circle fitter + ::LsqSums4 _szphi; // dphi / dz line fitter - XYZVectorF _center; - float _radius; + XYZVectorF _center; // circle fit results + float _radius; + float _circle_chisq_dof; // circle fit quality - // float _chi2; -//----------------------------------------------------------------------------- -// 2015-02-06 P.Murat: fit with non-equal weights - XY-only -//----------------------------------------------------------------------------- - // ::LsqSums4 _sxyw; - // XYZVectorF _cw; - // float _rw; - // float _chi2w; //----------------------------------------------------------------------------- // Z parameters; dfdz is the slope of phi vs z (=-sign(1.0,qBzdir)/(R*tandip)), // fz0 is the phi value of the particle where it goes through z=0 // note that dfdz has a physical ambiguity in q*zdir. //----------------------------------------------------------------------------- - float _dfdz; - float _fz0; + float _dfdz; // dphi / dz fit results + float _fz0; + float _dfdz_chisq_dof; // dphi / dz fit quality //----------------------------------------------------------------------------- // diagnostics, histogramming //----------------------------------------------------------------------------- Diag_t _diag; //----------------------------------------------------------------------------- -// structure used to organize thei strawHits for the pattern recognition +// structure used to organize the strawHits for the pattern recognition //----------------------------------------------------------------------------- // PanelZ_t _oTracker[kNTotalPanels]; // std::array _hitsUsed; diff --git a/CalPatRec/inc/CalHelixFinder_module.hh b/CalPatRec/inc/CalHelixFinder_module.hh index d9d04fd100..b2036e591d 100644 --- a/CalPatRec/inc/CalHelixFinder_module.hh +++ b/CalPatRec/inc/CalHelixFinder_module.hh @@ -75,7 +75,6 @@ namespace mu2e { // event object labels //----------------------------------------------------------------------------- std::string _shLabel ; // MakeStrawHit label (makeSH) - // std::string _shpLabel; std::string _timeclLabel; int _minNHitsTimeCluster; //min nhits within a TimeCluster after check of Delta-ray hits @@ -96,7 +95,6 @@ namespace mu2e { const ComboHitCollection* _chcol; const TimeClusterCollection* _timeclcol; - HelixTraj* _helTraj; CalHelixFinderAlg _hfinder; CalHelixFinderData _hfResult; std::vector _hels; // helicity values to fit @@ -161,6 +159,10 @@ namespace mu2e { int goodHitsTimeCluster(const TimeCluster* TimeCluster); void pickBestHelix(std::vector& HelVec, int &Index_best); + + void fillDiagnosticInfo(const std::vector& helix_seed_vec, + const std::vector& nHitsRatio_vec, + int index_best, const int nGoodTClusterHits); }; } #endif diff --git a/CalPatRec/inc/CalHelixFinder_types.hh b/CalPatRec/inc/CalHelixFinder_types.hh index 1338af1436..59141946d2 100644 --- a/CalPatRec/inc/CalHelixFinder_types.hh +++ b/CalPatRec/inc/CalHelixFinder_types.hh @@ -39,7 +39,6 @@ namespace mu2e { struct Data_t { const art::Event* event; - std::string shLabel; enum { kMaxSeeds = 100, kMaxHits = 200 }; diff --git a/CalPatRec/src/CalHelixFinderAlg.cc b/CalPatRec/src/CalHelixFinderAlg.cc index edaa872271..b0c8d3400e 100644 --- a/CalPatRec/src/CalHelixFinderAlg.cc +++ b/CalPatRec/src/CalHelixFinderAlg.cc @@ -318,11 +318,6 @@ namespace mu2e { //----------------------------------------------------------------------------- if (_diag > 0) saveResults(Helix, 0); - // if (_filter) { - // filterDist(Helix); - // if (_diag > 0) saveResults(_xyzp,Helix,1); - // } - doPatternRecognition(Helix); //--------------------------------------------------------------------------- // 2014-11-11 gianipez changed the following if() statement to test the @@ -331,7 +326,7 @@ namespace mu2e { //--------------------------------------------------------------------------- if (_debug != 0) { printf("[CalHelixFinderAlg::findHelix] Helix._nXYSh = %i goodPointsTrkCandidate = %i\n", - Helix._nXYSh, Helix._nStrawHits);//_goodPointsTrkCandidate); + Helix._nXYSh, Helix._nStrawHits); } if (Helix._nStrawHits < _minNHits ) { @@ -340,10 +335,12 @@ namespace mu2e { else if ((Helix._radius < _rmin) || (Helix._radius > _rmax)) { Helix._fit = TrkErrCode(TrkErrCode::fail,2); // initialization failure } - else if ((Helix._nXYSh < _minNHits) || (Helix._sxy.chi2DofCircle() > _chi2xyMax)) { + else if ((Helix._nXYSh < _minNHits) || (Helix._circle_chisq_dof > _chi2xyMax)) { Helix._fit = TrkErrCode(TrkErrCode::fail,3); // xy reconstruction failure } - else if ((Helix._nZPhiSh < _minNHits) || (Helix._szphi.chi2DofLine() > _chi2zphiMax) || + else if ((Helix._nZPhiSh < _minNHits) || + (Helix._dfdz_chisq_dof < 0.f) || // no valid line fit was found + (Helix._dfdz_chisq_dof > _chi2zphiMax) || (fabs(Helix._dfdz) < _minDfDz) || (fabs(Helix._dfdz) > _maxDfDz)) { Helix._fit = TrkErrCode(TrkErrCode::fail,4); // phi-z reconstruction failure } @@ -534,7 +531,7 @@ namespace mu2e { //----------------------------------------------------------------------------- for (int n=nmin; n<=nmax; n++) { const float x = dphidz + n*dx; - const int bin = std::min(nbinsX-1, int((x-minX)/stepX)); + const int bin = std::max(0, std::min(nbinsX-1, int((x-minX)/stepX))); hist[bin] += weight; } } @@ -1005,6 +1002,8 @@ namespace mu2e { } CHECK_RESIDUALS:; + constexpr bool endAtThreshold(false); // Whether or not to keep iterating down to the threshold or stop once it's acceptable + if (endAtThreshold && Helix._szphi.chi2DofLine() <= _chi2zphiMax) goto NEXT_ITERATION_END; dphi_max = _maxXDPhi; //reset the coordinates of the worst hit iworst.face = -1; @@ -1046,6 +1045,7 @@ namespace mu2e { goto NEXT_ITERATION; } } + NEXT_ITERATION_END:; } } //----------------------------------------------------------------------------- @@ -1063,6 +1063,7 @@ namespace mu2e { else if (success) { // update helix results Helix._fz0 = Helix._szphi.phi0(); Helix._dfdz = Helix._szphi.dfdz(); + Helix._dfdz_chisq_dof = Helix._szphi.chi2DofLine(); } if ((SeedIndex.face == 0 ) && (SeedIndex.panel == 0) && (SeedIndex.panelHitIndex == 0) && (_diag > 0)) { @@ -1151,6 +1152,9 @@ namespace mu2e { if (!success) { Helix._hitsUsed = hitsUsed; + // _szphi was cleared and partially rebuilt during the failed fit; + // invalidate the cached chi2 field so it cannot disagree with _szphi + Helix._dfdz_chisq_dof = -1.f; } return success; @@ -1461,6 +1465,7 @@ namespace mu2e { //----------------------------------------------------------------------------- Helix._center.SetXYZ(x0, y0, 0.0); Helix._radius = radius; + Helix._circle_chisq_dof = Helix._sxy.chi2DofCircle(); Helix._nComboHits += rescuedPoints; Helix._nStrawHits += rescuedStrawHits; @@ -1518,13 +1523,14 @@ namespace mu2e { void CalHelixFinderAlg::filterUsingPatternRecognition(CalHelixFinderData& Helix) { if (Helix._seedIndex.panel < 0) return; + if (Helix._dfdz == 0.) return; // undefined helix results int nActive(0), nActive_hel(0); int nSh = Helix._nFiltStrawHits; float straw_mean_radius(0), chi2_global_helix(0), total_weight(0); float x_center(Helix._center.x()), y_center(Helix._center.y()), radius(Helix._radius); - float fz0(Helix._fz0), lambda(Helix._dfdz != 0. ? 1./Helix._dfdz : 0.); + float fz0(Helix._fz0), dfdz(Helix._dfdz); XYZVectorF hel_pred(0., 0., 0.); PanelZ_t* panelz(0); @@ -1563,7 +1569,7 @@ namespace mu2e { float x = hit->pos().x(); float y = hit->pos().y(); float z = Helix._zFace[f]; - float phi_pred = fz0 + z/lambda; + float phi_pred = fz0 + z*dfdz; float x_pred = x_center + radius*cos(phi_pred); float y_pred = y_center + radius*sin(phi_pred); hel_pred.SetX(x_pred); @@ -1767,6 +1773,9 @@ namespace mu2e { rescueHits(Helix, HitInfo_t(StrawId::_ntotalfaces-1,0,-1), usePhiResid); if ((Helix._nXYSh - 1) != Helix._nZPhiSh) rc = doLinearFitPhiZ(Helix, HitInfo_t(StrawId::_ntotalfaces-1,0,-1), useIntelligentWeight);//the factor "-1" takes into account that the XY fit includes the target center + // rescueHits may have added points to _szphi directly (UsePhiResiduals==1); + // if doLinearFitPhiZ was skipped above, sync the field from the live fitter + else if(usePhiResid) Helix._dfdz_chisq_dof = Helix._szphi.chi2DofLine(); if (_debug != 0) printInfo(Helix); //-------------------------------------------------------------------------------------------------------------- @@ -1781,6 +1790,7 @@ namespace mu2e { if (rs == 1) { // update Helix Z-phi part Helix._dfdz = _hdfdz; Helix._fz0 = _hphi0; + Helix._dfdz_chisq_dof = -1.f; } } @@ -1790,6 +1800,9 @@ namespace mu2e { usePhiResid = 1; rescueHits(Helix, HitInfo_t(StrawId::_ntotalfaces-1,0,-1), usePhiResid); if ((Helix._nXYSh - 1) != Helix._nZPhiSh) rc = doLinearFitPhiZ(Helix, HitInfo_t(StrawId::_ntotalfaces-1,0,-1), useIntelligentWeight); //the factor "-1" takes into account that the XY fit includes the target center + // rescueHits may have added points to _szphi directly (UsePhiResiduals==1); + // if doLinearFitPhiZ was skipped above, sync the field from the live fitter + else if(usePhiResid) Helix._dfdz_chisq_dof = Helix._szphi.chi2DofLine(); if (_debug != 0) printInfo(Helix); strcpy(banner,"refineHelixParameters-after-doLinearFitPhiZ"); @@ -1921,9 +1934,12 @@ namespace mu2e { Helix._nXYSh = 0; Helix._nComboHits = 0; - Helix._sxy.addPoint(fCaloX,fCaloY,1./100.); - Helix._nXYSh += 1; - Helix._nComboHits += 1; + // Only add the calo cluster position if one is associated with this time cluster + if(fCaloTime > 0.) { + Helix._sxy.addPoint(fCaloX,fCaloY,1./100.); + Helix._nXYSh += 1; + Helix._nComboHits += 1; + } //------------------------------------------------------------------------------- // add stopping target center with a position error of 100 mm/sqrt(12) ~ 30mm => wt = 1/900 //------------------------------------------------------------------------------- @@ -1987,6 +2003,7 @@ namespace mu2e { Helix._center.SetX(Helix._sxy.x0()); Helix._center.SetY(Helix._sxy.y0()); Helix._radius = Helix._sxy.radius(); + Helix._circle_chisq_dof = Helix._sxy.chi2DofCircle(); if (_debug > 5) { @@ -2103,7 +2120,7 @@ namespace mu2e { chi2 = sxy.chi2DofCircle(); - if ((chi2 < chi2_min) || ( (i == SeedIndex.panelHitIndex) && (p == SeedIndex.panel) && (f == SeedIndex.face)) ) { + if (chi2 < chi2_min) { chi2_min = chi2; IWorst.face = f; IWorst.panel = p; @@ -2287,11 +2304,10 @@ namespace mu2e { //----------------------------------------------------------------------------- // update circle parameters //----------------------------------------------------------------------------- - // Trk._sxyw = sxyw; Trk._center.SetX(Trk._sxy.x0()); Trk._center.SetY(Trk._sxy.y0()); Trk._radius = Trk._sxy.radius(); - // Trk._chi2 = Trk._sxy.chi2DofCircle(); + Trk._circle_chisq_dof = Trk._sxy.chi2DofCircle(); }else { Trk._hitsUsed = hitsUsed; //restore the info of the used-hits that was originally passed to the procedure @@ -2590,7 +2606,7 @@ namespace mu2e { Helix._center.SetX(Helix._sxy.x0()); Helix._center.SetY(Helix._sxy.y0()); Helix._radius = Helix._sxy.radius(); - // Helix._chi2 = Helix._sxy.chi2DofCircle(); + Helix._circle_chisq_dof = Helix._sxy.chi2DofCircle(); F_END:; if (_debug > 5 ) { @@ -2945,10 +2961,12 @@ namespace mu2e { Helix._nXYSh = NPoints; Helix._radius = sxy.radius(); Helix._center.SetXYZ(sxy.x0(), sxy.y0(), 0.0); + Helix._circle_chisq_dof = sxy.chi2DofCircle(); Helix._nStrawHits = NPoints; Helix._nComboHits = NComboHits; Helix._dfdz = dfdz; Helix._fz0 = phi0 - dfdz*z_phi0; // *FLOAT_CHECK* + Helix._dfdz_chisq_dof = -1.f; radius_end = Helix._radius; //breakpoint -- radius before refine $ @@ -2958,11 +2976,8 @@ namespace mu2e { //----------------------------------------------------------------------------- if (rc < 0) return; //breakpoint -- radius after refine $ - // Helix._center.SetXYZ(Helix._cw.x(), Helix._cw.y(), 0.0); - // Helix._radius = Helix._rw; radius_end = Helix._radius; - // Helix._sxy = Helix._sxyw; // doWeightedCircleFit still adds the ST and the cluster Chi2 = Helix._sxy.chi2DofCircle(); NPoints = Helix._nStrawHits; // *FIXME* in principle, the fit can remove ST as well as the cluster @@ -2972,9 +2987,10 @@ namespace mu2e { // 2015-01-22 G. Pezzullo and P. Murat; update the dfdz value using all hits //----------------------------------------------------------------------------- int rs = findDfDz(Helix, SeedIndex); - if (rs ==1 ) { + if (rs ==1 ) { // 1 = success Helix._dfdz = _hdfdz; Helix._fz0 = _hphi0; + Helix._dfdz_chisq_dof = -1.f; // fill diag vector dfdzRes[1] = _hdfdz; dphi0Res[1] = _hphi0; @@ -2982,10 +2998,18 @@ namespace mu2e { //----------------------------------------------------------------------------- // 2015-01-23 G. Pezzu and P. Murat: when it fails, doLinearFitPhiZ returns negative value // in this case, use the previous value for dfdz and phi0 +// NOTE: _hphi0 is initialised to -9999. in resetTrackParamters() and is only updated by +// findDfDz() when >=2 stations have hits. If findDfDz() also failed, the fallback +// sets phi0_end = -9999., which is a sentinel value and NOT a valid phi0. +// The helix candidate is still stored in this state and compared against other +// candidates in searchBestTriplet; the final quality cuts in fitHelix will reject it +// only if the sentinel survives as _fz0 through to that check. +// A future improvement would be to return early from findTrack when rcPhiZ is false +// AND _hphi0 is still the sentinel value. //----------------------------------------------------------------------------- bool rcPhiZ = doLinearFitPhiZ(Helix, SeedIndex); - if (rcPhiZ) { + if (rcPhiZ) { // fit was successful dfdz_end = Helix._dfdz; phi0_end = Helix._fz0; // fill diagnostic vector @@ -2998,27 +3022,6 @@ namespace mu2e { //FIXME! implement a function countUsedHits(Helix, SeedIndex, NComboHits, NPoints); - // for (int f=SeedIndex.face; fpanelZs[p]; - // int nhitsPerPanel = panelz->nChHits(); - // int seedPanelIndex(0); - // if (nhitsPerPanel == 0) continue; - // if ( (f==SeedIndex.face) && (p==SeedIndex.panel) && (SeedIndex.panelHitIndex >=0)) seedPanelIndex = SeedIndex.panelHitIndex - panelz->idChBegin; - - // for (int i=seedPanelIndex; iidChBegin + i; - // hit = &Helix._chHitsToProcess[index]; - // if (Helix._hitsUsed[index] > 0 ) { - // ++NComboHits; - // NPoints += hit->nStrawHits(); - // } - // } - // }//endl panels loop - // } } else { dfdz_end = _hdfdz; @@ -3075,8 +3078,12 @@ namespace mu2e { //----------------------------------------------------------------------------- Helix._center.SetXYZ(Helix._sxy.x0(),Helix._sxy.y0(), 0.0); Helix._radius = radius_end; + Helix._circle_chisq_dof = Helix._sxy.chi2DofCircle(); Helix._fz0 = phi0_end; Helix._dfdz = dfdz_end; + // _dfdz_chisq_dof is already set by doLinearFitPhiZ when rcPhiZ==true; + // only reset it to 0 when the fit failed (histogram estimate, no valid chi2) + if (!rcPhiZ) Helix._dfdz_chisq_dof = -1.f; if (_diag > 0){ Helix._diag.loopId_4 = _findTrackLoopIndex; @@ -3106,7 +3113,7 @@ namespace mu2e { int firstPanel(0); if (f == SeedIndex.face) firstPanel = SeedIndex.panel; for (int p=firstPanel; ppanelZs[p];//Helix._oTracker[p]; + panelz = &facez->panelZs[p]; int nhitsPerPanel = panelz->nChHits(); int seedPanelIndex(0); if (nhitsPerPanel == 0) continue; @@ -3116,17 +3123,8 @@ namespace mu2e { int index = panelz->idChBegin + i; hit = &Helix._chHitsToProcess[index]; - // int index = facez->evalUniqueHitIndex(f,p,i);//p*FaceZ_t::kNMaxHitsPerPanel + i; if (Helix._hitsUsed[index] != 1) continue; - // if (j < Helix.maxIndex()) { - // Helix._diag.dist[j] = hit->_drFromPred; - // Helix._diag.dz [j] = hit->_dzFromSeed; - // ++j; - // } - // else { - // printf("ERROR in CalHelixFinderAlg::findTrack : index out limits. IGNORE; \n"); - // } }//end loop over the hits within a panel }//end panel loop }//end face loop diff --git a/CalPatRec/src/CalHelixFinderData.cc b/CalPatRec/src/CalHelixFinderData.cc index 7b91201e47..d2126ffa35 100644 --- a/CalPatRec/src/CalHelixFinderData.cc +++ b/CalPatRec/src/CalHelixFinderData.cc @@ -113,9 +113,11 @@ namespace mu2e { _szphi.clear(); _radius = -1.; + _circle_chisq_dof = 1e10; _dfdz = -1.e6; _fz0 = -1.e6; + _dfdz_chisq_dof = 1e10; _nFiltPoints = 0; _nFiltStrawHits = 0; @@ -185,9 +187,11 @@ void CalHelixFinderData::clearHelixInfo() { _szphi.clear(); _radius = -1.; + _circle_chisq_dof = 1e10; _dfdz = -1.e6; _fz0 = -1.e6; + _dfdz_chisq_dof = 1e10; _nXYSh = 0; _nZPhiSh = 0; @@ -216,9 +220,11 @@ void CalHelixFinderData::clearHelixInfo() { _szphi.clear(); // _chi2 = -1.; _radius = -1.; + _circle_chisq_dof = 1e10; _dfdz = -1.e6; _fz0 = -1.e6; + _dfdz_chisq_dof = 1e10; _nXYSh = 0; diff --git a/CalPatRec/src/CalHelixFinder_module.cc b/CalPatRec/src/CalHelixFinder_module.cc index b7ec9a973b..ecad489323 100644 --- a/CalPatRec/src/CalHelixFinder_module.cc +++ b/CalPatRec/src/CalHelixFinder_module.cc @@ -20,7 +20,6 @@ #include "Offline/GeometryService/inc/DetectorSystem.hh" #include "Offline/CalorimeterGeom/inc/DiskCalorimeter.hh" #include "Offline/ConfigTools/inc/ConfigFileLookupPolicy.hh" -// #include "CalPatRec/inc/KalFitResult.hh" #include "Offline/RecoDataProducts/inc/StrawHitIndex.hh" #include "Offline/RecoDataProducts/inc/HelixHit.hh" @@ -41,6 +40,8 @@ #include "TSystem.h" #include "TInterpreter.h" +#include + using namespace boost::accumulators; using CLHEP::HepVector; using CLHEP::Hep3Vector; @@ -65,11 +66,14 @@ namespace mu2e { consumes(_shLabel); consumes(_timeclLabel); + // Initialize the list of helicities to consider std::vector helvals = config().Helicities(); for(auto hv : helvals) { - Helicity hel(hv); - _hels.push_back(hel); + Helicity hel(hv); + _hels.push_back(hel); } + + // In single output, put all helices by all helicities into one helix collection, otherwise produce a collection per helicity if (_doSingleOutput){ produces(); } else { @@ -79,12 +83,9 @@ namespace mu2e { } _tpart = (TrkParticle::type)_fitparticle; - //----------------------------------------------------------------------------- - // provide for interactive disanostics - //----------------------------------------------------------------------------- - _helTraj = 0; - - _data.shLabel = _shLabel; +//----------------------------------------------------------------------------- +// provide for interactive diagnostics +//----------------------------------------------------------------------------- if (_debugLevel != 0) _printfreq = 1; @@ -92,20 +93,21 @@ namespace mu2e { else _hmanager = std::make_unique(); } - //----------------------------------------------------------------------------- - // destructor - //----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +// destructor +//----------------------------------------------------------------------------- CalHelixFinder::~CalHelixFinder() { - if (_helTraj) delete _helTraj; } - //----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- void CalHelixFinder::beginJob(){ - art::ServiceHandle tfs; - _hmanager->bookHistograms(tfs); + if(_diagLevel != 0) { // only book histograms if running diagnostics + art::ServiceHandle tfs; + _hmanager->bookHistograms(tfs); + } } - //----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- void CalHelixFinder::beginRun(art::Run& ) { mu2e::GeomHandle bfmgr; mu2e::GeomHandle det; @@ -168,9 +170,9 @@ namespace mu2e { } - //----------------------------------------------------------------------------- - // find the input data objects - //----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +// find the input data objects +//----------------------------------------------------------------------------- bool CalHelixFinder::findData(const art::Event& evt) { if (evt.getByLabel(_shLabel, _strawhitsH)) { @@ -179,41 +181,30 @@ namespace mu2e { else { _chcol = nullptr; printf(" >>> ERROR in CalHelixFinder::findData: StrawHitCollection with label=%s not found.\n", - _shLabel.data()); + _shLabel.data()); } - // art::Handle shposH; - // if (evt.getByLabel(_shpLabel,shposH)) { - // _shpcol = shposH.product(); - // } - // else { - // _shpcol = 0; - // printf(" >>> ERROR in CalHelixFinder::findData: StrawHitPositionCollection with label=%s not found.\n", - // _shpLabel.data()); - // } - if (evt.getByLabel(_timeclLabel, _timeclcolH)) { _timeclcol = _timeclcolH.product(); } else { _timeclcol = nullptr; printf(" >>> ERROR in CalHelixFinder::findData: TimeClusterCollection with label=%s not found.\n", - _timeclLabel.data()); + _timeclLabel.data()); } - //----------------------------------------------------------------------------- - // done - //----------------------------------------------------------------------------- - return (_chcol != nullptr) //&& (_shpcol != nullptr) - && (_timeclcol != nullptr); +//----------------------------------------------------------------------------- +// done +//----------------------------------------------------------------------------- + return (_chcol != nullptr) && (_timeclcol != nullptr); } - //----------------------------------------------------------------------------- - // event entry point - //----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +// event entry point +//----------------------------------------------------------------------------- void CalHelixFinder::produce(art::Event& event ) { const char* oname = "CalHelixFinder::filter"; - // diagnostic info + // diagnostic info _data.event = &event; _iev = event.id().event(); int nGoodTClusterHits(0); @@ -228,26 +219,24 @@ namespace mu2e { _data.nseeds [counter] = 0; ++counter; } - }else { + } else { helcols[0] = std::unique_ptr(new HelixSeedCollection()); _data.nseeds [counter] = 0; } - // unique_ptr outseeds(new HelixSeedCollection); - //----------------------------------------------------------------------------- - // find the data - //----------------------------------------------------------------------------- + +//----------------------------------------------------------------------------- +// find the data +//----------------------------------------------------------------------------- if (!findData(event)) { printf("%s ERROR: No straw hits found, RETURN\n", oname); - goto END; + goto END; } - //----------------------------------------------------------------------------- - // loop over found time peaks - for us, - "eligible" calorimeter clusters - //----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +// loop over found time peaks - for us, - "eligible" calorimeter clusters +//----------------------------------------------------------------------------- _hfResult._tpart = _tpart; _hfResult._fdir = _fdir; _hfResult._chcol = _chcol; - // _hfResult._shpos = _shpcol; - //_hfResult._shfcol = _shfcol; _data.nTimePeaks = _timeclcol->size(); for (int ipeak=0; ipeak<_data.nTimePeaks; ipeak++) { @@ -258,28 +247,28 @@ namespace mu2e { std::vector helix_seed_vec; - //----------------------------------------------------------------------------- - // create track definitions for the helix fit from this initial information - // track fitting objects for this peak - //----------------------------------------------------------------------------- - _hfResult.clearTempVariables();//clearTimeClusterInfo(); +//----------------------------------------------------------------------------- +// create track definitions for the helix fit from this initial information +// track fitting objects for this peak +//----------------------------------------------------------------------------- + _hfResult.clearTempVariables(); _hfResult._timeCluster = tc; _hfResult._timeClusterPtr = art::Ptr(_timeclcolH,ipeak); - //----------------------------------------------------------------------------- - // fill the face-order hits collector - //----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +// fill the face-order hits collector +//----------------------------------------------------------------------------- _hfinder.fillFaceOrderedHits(_hfResult); - //----------------------------------------------------------------------------- - // Step 1: now loop over the two possible helicities. - // Find initial helical approximation of a track for both hypothesis - //----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +// Step 1: now loop over the two possible helicities. +// Find initial helical approximation of a track for both hypothesis +//----------------------------------------------------------------------------- std::vector nHitsRatio_vec; for (size_t i=0; i<_hels.size(); ++i){ - //----------------------------------------------------------------------------- - // create track definitions for the helix fit from this initial information - // track fitting objects for this peak - //----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +// create track definitions for the helix fit from this initial information +// track fitting objects for this peak +//----------------------------------------------------------------------------- CalHelixFinderData tmpResult(_hfResult); tmpResult.clearHelixInfo(); @@ -290,7 +279,7 @@ namespace mu2e { if (!rc) continue; if(!tmpResult.helix()) { std::cout << "[CalHelixFinder::" << __func__ << "] " << event.id() - << " Helix found but nullptr returned!!\n"; + << " Helix found but nullptr returned!!\n"; continue; } HelixSeed tmp_helix_seed; @@ -305,15 +294,15 @@ namespace mu2e { if (helix_seed_vec.size() == 0) continue; - //----------------------------------------------------------------------------- - // now select the best helix to avoid duplicates - //----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +// now select the best helix to avoid duplicates +//----------------------------------------------------------------------------- int index_best(-1); pickBestHelix(helix_seed_vec, index_best); - //----------------------------------------------------------------------------- - // fill seed information - //----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +// fill seed information +//----------------------------------------------------------------------------- if ( (index_best>=0) && (index_best < 2) ){ Helicity hel_best = helix_seed_vec[index_best]._helix._helicity; if (_doSingleOutput) { @@ -335,106 +324,20 @@ namespace mu2e { } } - // helix_seed_vec[index_best]._status.merge(TrkFitFlag::helixOK); - // outseeds->push_back(helix_seed_vec[index_best]); - if (_diagLevel > 0) { - //-------------------------------------------------------------------------------- - // fill diagnostic information - //-------------------------------------------------------------------------------- - int nhitsMin(15); - double mm2MeV = (3/10.)*_bz0; - - int loc = _data.nseeds[0]; - if (loc < _data.maxSeeds()) { - if (index_best == 2){ - index_best = 0; - } - int nhits = helix_seed_vec[index_best]._hhits.size(); - _data.ntclhits[loc]= nGoodTClusterHits; - _data.nhits[loc] = nhits; - _data.radius[loc] = helix_seed_vec[index_best].helix().radius(); - _data.pT[loc] = mm2MeV*_data.radius[loc]; - _data.p[loc] = _data.pT[loc]/std::cos( std::atan(helix_seed_vec[index_best].helix().lambda()/_data.radius[loc])); - - _data.chi2XY[loc] = _hfResult._sxy.chi2DofCircle(); - _data.chi2ZPhi[loc] = _hfResult._szphi.chi2DofLine(); - - _data.nHitsRatio[loc] = nHitsRatio_vec[index_best]; - _data.eDepAvg[loc] = helix_seed_vec[index_best]._eDepAvg; - - _data.nseeds[0]++; - _data.good[loc] = 0; - if (nhits >= nhitsMin) { - _data.nseeds[1]++; - _data.good[loc] = 1; - } - _data.nStationPairs[loc] = _hfResult._diag.nStationPairs; - - _data.dr [loc] = _hfResult._diag.dr; - _data.shmeanr [loc] = _hfResult._diag.straw_mean_radius; - _data.chi2d_helix [loc] = _hfResult._diag.chi2d_helix; - if (_hfResult._diag.chi2d_helix>3) printf("[%s] : chi2Helix = %10.3f event number %8i\n", oname,_hfResult._diag.chi2d_helix,_iev); - //----------------------------------------------------------------------------- - // info of the track candidate after the first loop with findtrack on CalHelixFinderAlg::doPatternRecognition - //----------------------------------------------------------------------------- - _data.loopId [loc] = _hfResult._diag.loopId_4; - if (_hfResult._diag.loopId_4 == 1) { - _data.chi2d_loop0 [loc] = _hfResult._diag.chi2_dof_circle_12; - _data.chi2d_line_loop0 [loc] = _hfResult._diag.chi2_dof_line_13; - _data.npoints_loop0 [loc] = _hfResult._diag.n_active_11; - - } - if (_hfResult._diag.loopId_4 == 2){ - _data.chi2d_loop1 [loc] = _hfResult._diag.chi2_dof_circle_12; - _data.chi2d_line_loop1 [loc] = _hfResult._diag.chi2_dof_line_13; - _data.npoints_loop1 [loc] = _hfResult._diag.n_active_11; - } - - //-------------------------------------------------------------------------------- - // info of the track candidate during the CAlHelixFinderAlg::findTrack loop - //-------------------------------------------------------------------------------- - int counter(0); - for (unsigned i=0; i<_hfResult._hitsUsed.size(); ++i){ - if (_hfResult._hitsUsed[i] != 1) continue; - ++counter; - } - // for (int f=0; fpanelZs[p];//&_hfResult._oTracker[p]; - // int nhits = panelz->fNHits; - // if (nhits == 0) continue; - - // for (int i=0; i_chHitsToProcess.at(i); - // int index = facez->evalUniqueHitIndex(f,p,i);//p*CalHelixFinderData::kNMaxHitsPerPanel + i; - // if (_hfResult._hitsUsed[index] != 1) continue; - - // // double dzFromSeed = hit->_dzFromSeed; //distance form the hit used to seed the 3D-search - // // double drFromPred = hit->_drFromPred; //distance from prediction - // // _data.hitDzSeed[loc][counter] = dzFromSeed; - // // _data.hitDrPred[loc][counter] = drFromPred; - // ++counter; - // }//end loop over the hits within a panel - // }//end panels loop - // }//end faces loop - } - else { - printf(" N(seeds) > %i, IGNORE SEED\n",_data.maxSeeds()); - } - } - - } - //-------------------------------------------------------------------------------- - // fill histograms - //-------------------------------------------------------------------------------- if (_diagLevel > 0) { - _hmanager->fillHistograms(&_data); + fillDiagnosticInfo(helix_seed_vec, nHitsRatio_vec, + index_best, nGoodTClusterHits); } - //----------------------------------------------------------------------------- - // put reconstructed tracks into the event record - //----------------------------------------------------------------------------- -END:; + } +//-------------------------------------------------------------------------------- +// fill histograms +//-------------------------------------------------------------------------------- + if (_diagLevel > 0) _hmanager->fillHistograms(&_data); + +//----------------------------------------------------------------------------- +// put reconstructed tracks into the event record +//----------------------------------------------------------------------------- + END:; int nseeds(0); if (_doSingleOutput) { nseeds += helcols[0]->size(); @@ -445,258 +348,290 @@ END:; event.put(std::move(helcols[hel]),Helicity::name(hel)); } } - } + } - //----------------------------------------------------------------------------- - // - //----------------------------------------------------------------------------- - void CalHelixFinder::endJob() { - // does this cause the file to close? - art::ServiceHandle tfs; - } +//----------------------------------------------------------------------------- +// +//----------------------------------------------------------------------------- + void CalHelixFinder::endJob() { + } - //-------------------------------------------------------------------------------- - // set helix parameters - //----------------------------------------------------------------------------- - void CalHelixFinder::initHelixSeed(HelixSeed& HelSeed, CalHelixFinderData& HfResult) { +//-------------------------------------------------------------------------------- +// set helix parameters +//----------------------------------------------------------------------------- + void CalHelixFinder::initHelixSeed(HelixSeed& HelSeed, CalHelixFinderData& HfResult) { - HelixTraj* hel = HfResult.helix(); + HelixTraj* hel = HfResult.helix(); - double helixRadius = 1./fabs(hel->omega()); - double impactParam = hel->d0(); - double phi0 = hel->phi0(); - // double x0 = -(helixRadius + impactParam)*sin(phi0)*sig; - // double y0 = (helixRadius + impactParam)*cos(phi0)*sig; + double helixRadius = 1./fabs(hel->omega()); + double impactParam = hel->d0(); + double phi0 = hel->phi0(); - double x0 = -(1/hel->omega()+impactParam)*sin(phi0); - double y0 = (1/hel->omega()+impactParam)*cos(phi0); + double x0 = -(1/hel->omega()+impactParam)*sin(phi0); + double y0 = (1/hel->omega()+impactParam)*cos(phi0); - double dfdz = 1./hel->tanDip()/helixRadius; - double z0 = hel->z0(); - // center of the helix in the transverse plane - Hep3Vector center(x0, y0, 0); - //define the reconstructed helix parameters - HelSeed._helix._rcent = center.perp(); - HelSeed._helix._fcent = center.phi(); - HelSeed._helix._radius = helixRadius; - HelSeed._helix._lambda = 1./dfdz*_hfinder._dfdzsign; + double dfdz = 1./hel->tanDip()/helixRadius; + double z0 = hel->z0(); + // center of the helix in the transverse plane + Hep3Vector center(x0, y0, 0); + //define the reconstructed helix parameters + HelSeed._helix._rcent = center.perp(); + HelSeed._helix._fcent = center.phi(); + HelSeed._helix._radius = helixRadius; + HelSeed._helix._lambda = 1./dfdz*_hfinder._dfdzsign; - HelSeed._helix._fz0 = phi0 - M_PI/2.*_hfinder._dfdzsign -z0*hel->omega()/hel->tanDip() ; + HelSeed._helix._fz0 = phi0 - M_PI/2.*_hfinder._dfdzsign -z0*hel->omega()/hel->tanDip() ; - HelSeed._helix._helicity = HfResult._helicity; - HelSeed._status.merge(TrkFitFlag::CPRHelix); + HelSeed._helix._helicity = HfResult._helicity; + HelSeed._status.merge(TrkFitFlag::CPRHelix); - //include also the values of the chi2d - HelSeed._helix._chi2dXY = HfResult._sxy.chi2DofCircle(); - HelSeed._helix._chi2dZPhi = HfResult._szphi.chi2DofLine(); + //include also the values of the chi2d + HelSeed._helix._chi2dXY = HfResult._sxy.chi2DofCircle(); + HelSeed._helix._chi2dZPhi = HfResult._szphi.chi2DofLine(); - //now evaluate the helix T0 using the calorimeter cluster - double mm2MeV = (3/10.)*_bz0; - double tandip = hel->tanDip(); - double mom = helixRadius*mm2MeV/std::cos( std::atan(tandip)); - double beta = _tpart.beta(mom); + //now evaluate the helix T0 using the calorimeter cluster + double mm2MeV = (3/10.)*_bz0; + double tandip = hel->tanDip(); + double mom = helixRadius*mm2MeV/std::cos( std::atan(tandip)); + double beta = _tpart.beta(mom); - double hel_t0; - double pitchAngle = M_PI/2. - std::atan(tandip); + double hel_t0; + double pitchAngle = M_PI/2. - std::atan(tandip); - HelSeed._timeCluster = HfResult._timeClusterPtr; + HelSeed._timeCluster = HfResult._timeClusterPtr; - // cluster hits assigned to the reconsturcted Helix + // cluster hits assigned to the reconsturcted Helix - int nhits = HfResult.nGoodHits(); - // printf("[CalHelixFinder::initHelixSeed] radius = %2.3f x0 = %2.3f y0 = %2.3f dfdz = %2.3e nhits = %i chi2XY = %2.3f chi2PHIZ = %2.3f\n", - // helixRadius, center.x(), center.y(), dfdz, nhits, HfResult._sxyw.chi2DofCircle(), HfResult._srphi.chi2DofLine()); - // printf("[CalHelixFinder::initHelixSeed] Index X Y Z PHI\n"); + int nhits = HfResult.nGoodHits(); - // double z_start(0); - HelSeed._hhits.setParent(_chcol->parent()); + HelSeed._hhits.setParent(_chcol->parent()); - for (int i=0; i_chHitsToProcess.at(hitInfo->panelHitIndex); - ComboHit hhit(*hit); - HelSeed._hhits.push_back(hhit); - } + for (int i=0; i_chHitsToProcess.at(hitInfo->panelHitIndex); + ComboHit hhit(*hit); + HelSeed._hhits.push_back(hhit); + } + + if(HfResult._timeClusterPtr->hasCaloCluster()) { + CLHEP::Hep3Vector gpos,tpos; + + gpos = _hfinder._calorimeter->geomUtil().diskToMu2e(HfResult._timeClusterPtr->caloCluster()->diskID(), + HfResult._timeClusterPtr->caloCluster()->cog3Vector()); + tpos = _hfinder._calorimeter->geomUtil().mu2eToTracker(gpos); + + hel_t0 = HfResult._timeClusterPtr->caloCluster()->time() - (tpos.z() - z0)/sin(pitchAngle)/(beta*CLHEP::c_light); - if(HfResult._timeClusterPtr->hasCaloCluster()) { - CLHEP::Hep3Vector gpos,tpos; + } + else { - gpos = _hfinder._calorimeter->geomUtil().diskToMu2e(HfResult._timeClusterPtr->caloCluster()->diskID(), - HfResult._timeClusterPtr->caloCluster()->cog3Vector()); - tpos = _hfinder._calorimeter->geomUtil().mu2eToTracker(gpos); + double tSum = 0.; - hel_t0 = HfResult._timeClusterPtr->caloCluster()->time() - (tpos.z() - z0)/sin(pitchAngle)/(beta*CLHEP::c_light); + size_t nHits = HelSeed.hits().size(); + if(nHits == 0) { + hel_t0 = 0.; } else { - double tSum = 0.; + for(size_t i = 0; inhits(); - int ngoodhits(0); - // double minT(500.), maxT(2000.); - for (int i=0; ihits().at(i); - // StrawHitFlag flag = _shfcol->at(index); - ComboHit sh = _chcol ->at(index); - StrawHitFlag flag = sh.flag(); - int bkg_hit = flag.hasAnyProperty(StrawHitFlag::bkg); - if (bkg_hit) continue; - // if ( (sh.time() < minT) || (sh.time() > maxT) ) continue; - - ngoodhits += sh.nStrawHits(); - } +//----------------------------------------------------------------------------- + int CalHelixFinder::goodHitsTimeCluster(const TimeCluster* TCluster){ + int nhits = TCluster->nhits(); + int ngoodhits(0); + for (int i=0; ihits().at(i); + ComboHit sh = _chcol ->at(index); + StrawHitFlag flag = sh.flag(); + int bkg_hit = flag.hasAnyProperty(StrawHitFlag::bkg); + if (bkg_hit) continue; + + ngoodhits += sh.nStrawHits(); + } + + return ngoodhits; + } - return ngoodhits; +//-------------------------------------------------------------------------------- +// function to select the best Helix among the results of the two helicity hypo +//-------------------------------------------------------------------------------- + void CalHelixFinder::pickBestHelix(std::vector& HelVec, int &Index_best){ + if (HelVec.size() == 1) { + Index_best = 0; + return; } - //-------------------------------------------------------------------------------- - // function to select the best Helix among the results of the two helicity hypo - //-------------------------------------------------------------------------------- - void CalHelixFinder::pickBestHelix(std::vector& HelVec, int &Index_best){ - if (HelVec.size() == 1) { - Index_best = 0; - return; - } +//----------------------------------------------------------------------------- +// at Mu2e, 2 helices with different helicity could be duplicates of each other +//----------------------------------------------------------------------------- - const HelixSeed *h1, *h2; - const ComboHitCollection *tlist, *clist; - int nh1, nh2, natc(0); - const mu2e::HelixHit *hitt, *hitc; + const HelixSeed *h1, *h2; + const ComboHitCollection *tlist, *clist; + int nh1, nh2, natc(0); - h1 = &HelVec[0]; - //------------------------------------------------------------------------------ - // check if an AlgorithmID collection has been created by the process - //----------------------------------------------------------------------------- - tlist = &h1->hits(); - nh1 = tlist->size(); + constexpr float overlap_threshold(0.5f); // 50% overlap to match helices - h2 = &HelVec[1]; - //----------------------------------------------------------------------------- - // at Mu2e, 2 helices with different helicity could be duplicates of each other - //----------------------------------------------------------------------------- - clist = &h2->hits(); - nh2 = clist->size(); + // Helix 1 + h1 = &HelVec[0]; + tlist = &h1->hits(); + nh1 = tlist->size(); - //----------------------------------------------------------------------------- - // check the number of common hits - //----------------------------------------------------------------------------- - for (int k=0; kat(k); - for (int l=0; lat(l); - if (hitt->index() == hitc->index()) { - natc += 1; - break; - } - } - } + // Helix 2 + h2 = &HelVec[1]; + clist = &h2->hits(); + nh2 = clist->size(); +//----------------------------------------------------------------------------- +// check the number of common hits +//----------------------------------------------------------------------------- + std::unordered_set hits_1; // unordered set for O(1) lookups + hits_1.reserve(nh1); - if ((natc > nh1/2.) || (natc > nh2/2.)) { + // Fill the set + for(int k = 0; k < nh1; ++k) { + hits_1.insert(tlist->at(k).index()); + } - //----------------------------------------------------------------------------- - // pick the helix with the largest number of hits - //----------------------------------------------------------------------------- - if (nh2 > nh1) { - //----------------------------------------------------------------------------- - // h2 is a winner, no need to save h1 - //----------------------------------------------------------------------------- - Index_best = 1; - return; - } - else if (nh1 > nh2){ - //----------------------------------------------------------------------------- - // h1 is a winner, mark h2 in hope that it will be OK, continue looping - //----------------------------------------------------------------------------- - Index_best = 0; - return; - } - //----------------------------------------------------------------------------- - // in case they have the exact amount of hits, pick the one with better chi2dZphi - //----------------------------------------------------------------------------- - if (nh1 == nh2) { - float chi2dZphi_h1 = h1->helix().chi2dZPhi(); - float chi2dZphi_h2 = h2->helix().chi2dZPhi(); - if (chi2dZphi_h1 < chi2dZphi_h2){ - Index_best = 0; - return; - }else { - Index_best = 1; - return; - } - } - }else { - //----------------------------------------------------------------------------- - // this is the case where we consider the two helices independent, so we want - // to store both - //----------------------------------------------------------------------------- - Index_best = 2; - return; - } + // Count overlapping hits + for(int l = 0; l < nh2; l++) { + if(hits_1.contains(clist->at(l).index())) ++natc; + } + // Check if at least one shares half of its hits with the other helix + if ((natc > nh1*overlap_threshold) || (natc > nh2*overlap_threshold)) { + + //----------------------------------------------------------------------------- + // pick the helix with the largest number of hits + //----------------------------------------------------------------------------- + if (nh2 > nh1) Index_best = 1; // helix two has more hits + else if (nh1 > nh2) Index_best = 0; // helix one has more hits + else Index_best = (h1->helix().chi2dZPhi() < h2->helix().chi2dZPhi()) ? 0 : 1; // equal hits --> pick the best chi^2 + } else { +//----------------------------------------------------------------------------- +// this is the case where we consider the two helices independent, so we want +// to store both +//----------------------------------------------------------------------------- + Index_best = 2; } } - using mu2e::CalHelixFinder; - DEFINE_ART_MODULE(CalHelixFinder) + +//----------------------------------------------------------------------------- +// fill diagnostic information +//----------------------------------------------------------------------------- + void CalHelixFinder::fillDiagnosticInfo(const std::vector& helix_seed_vec, + const std::vector& nHitsRatio_vec, + int index_best, const int nGoodTClusterHits) { + const char* oname = "CalHelixFinder::fillDiagnosticsInfo"; + int nhitsMin(15); + double mm2MeV = (3/10.)*_bz0; + + int loc = _data.nseeds[0]; + if (loc < _data.maxSeeds()) { + if (index_best == 2){ + index_best = 0; + } + int nhits = helix_seed_vec[index_best]._hhits.size(); + _data.ntclhits[loc]= nGoodTClusterHits; + _data.nhits[loc] = nhits; + _data.radius[loc] = helix_seed_vec[index_best].helix().radius(); + _data.pT[loc] = mm2MeV*_data.radius[loc]; + _data.p[loc] = _data.pT[loc]/std::cos( std::atan(helix_seed_vec[index_best].helix().lambda()/_data.radius[loc])); + + _data.chi2XY[loc] = _hfResult._sxy.chi2DofCircle(); + _data.chi2ZPhi[loc] = _hfResult._szphi.chi2DofLine(); + + _data.nHitsRatio[loc] = nHitsRatio_vec[index_best]; + _data.eDepAvg[loc] = helix_seed_vec[index_best]._eDepAvg; + + _data.nseeds[0]++; + _data.good[loc] = 0; + if (nhits >= nhitsMin) { + _data.nseeds[1]++; + _data.good[loc] = 1; + } + _data.nStationPairs[loc] = _hfResult._diag.nStationPairs; + + _data.dr [loc] = _hfResult._diag.dr; + _data.shmeanr [loc] = _hfResult._diag.straw_mean_radius; + _data.chi2d_helix [loc] = _hfResult._diag.chi2d_helix; + if (_hfResult._diag.chi2d_helix>3) printf("[%s] : chi2Helix = %10.3f event number %8i\n", oname,_hfResult._diag.chi2d_helix,_iev); + //----------------------------------------------------------------------------- + // info of the track candidate after the first loop with findtrack on CalHelixFinderAlg::doPatternRecognition + //----------------------------------------------------------------------------- + _data.loopId [loc] = _hfResult._diag.loopId_4; + if (_hfResult._diag.loopId_4 == 1) { + _data.chi2d_loop0 [loc] = _hfResult._diag.chi2_dof_circle_12; + _data.chi2d_line_loop0 [loc] = _hfResult._diag.chi2_dof_line_13; + _data.npoints_loop0 [loc] = _hfResult._diag.n_active_11; + + } + if (_hfResult._diag.loopId_4 == 2){ + _data.chi2d_loop1 [loc] = _hfResult._diag.chi2_dof_circle_12; + _data.chi2d_line_loop1 [loc] = _hfResult._diag.chi2_dof_line_13; + _data.npoints_loop1 [loc] = _hfResult._diag.n_active_11; + } + + //-------------------------------------------------------------------------------- + // info of the track candidate during the CAlHelixFinderAlg::findTrack loop + //-------------------------------------------------------------------------------- + int counter(0); + for (unsigned i=0; i<_hfResult._hitsUsed.size(); ++i){ + if (_hfResult._hitsUsed[i] != 1) continue; + ++counter; + } + } + else { + printf(" N(seeds) > %i, IGNORE SEED\n",_data.maxSeeds()); + } + } +} + +using mu2e::CalHelixFinder; +DEFINE_ART_MODULE(CalHelixFinder) diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index cc430349fc..f46290275d 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -38,17 +38,14 @@ namespace mu2e { template bool extrapolateST(KKTrack& ktrk,TimeDir trkdir) const; template bool extrapolateTracker(KKTrack& ktrk,TimeDir tdir) const; template bool extrapolateTSDA(KKTrack& ktrk,TimeDir tdir) const; - template bool extrapolateTCRV(KKTrack& ktrk,TimeDir tdir) const; template void toOPA(KKTrack& ktrk, double tstart, TimeDir tdir) const; - template void toCRV(KKTrack& ktrk, TimeDir tdir) const; private: int debug_; - double btol_, intertol_, maxdt_, minv_; - bool backToTracker_, toOPA_, toTrackerEnds_, upstream_, toCRV_; + double btol_, intertol_, maxdt_; + bool backToTracker_, toOPA_, toTrackerEnds_, upstream_; double ipathick_ = 0.511; // ipa thickness: should come from geometry service TODO double stthick_ = 0.1056; // st foil thickness: should come from geometry service TODO - double tcrvthick_ = 150.0; // test CRV sector thickness: should come from geometry service TODO }; KKExtrap::KKExtrap(KKExtrapConfig const& extrapconfig) : @@ -56,10 +53,8 @@ namespace mu2e { btol_(extrapconfig.btol()), intertol_(extrapconfig.interTol()), maxdt_(extrapconfig.MaxDt()), - minv_(extrapconfig.MinV()), backToTracker_(extrapconfig.BackToTracker()), toOPA_(extrapconfig.ToOPA()), - toCRV_(extrapconfig.ToCRV()), toTrackerEnds_(extrapconfig.ToTrackerEnds()), upstream_(extrapconfig.Upstream()) {} @@ -100,8 +95,6 @@ namespace mu2e { } // optionally test for intersection with the OPA if(toOPA_)toOPA(ktrk,starttime,tdir); - // optionally test for intersection with the CRV - if(toCRV_)toCRV(ktrk,starttime,tdir); } } @@ -278,40 +271,5 @@ namespace mu2e { ktrk.addIntersection(OPASID,opainter); } } - - template void KKExtrap::toTCRV(KKTrack& ktrk,TimeDir tdir) const { - GeomHandle kkg_h; - GeomHandle kkmat_h; - // extrapolate to the extracted CRV - auto TCRV = ExtrapolateTCRV(maxdt_,btol_,intertol_,minv_,*kkg_h->TCRV(),debug_); - - auto const& ftraj = ktrk.fitTraj(); - static const SurfaceId TCRVSID("TCRV"); - auto dir0 = ftraj.direction(ftraj.t0()); - TimeDir tdir = (dir0.Y() > 0) ? TimeDir::forwards : TimeDir::backwards; - double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); - bool hadintersection = false; - do { - // iterate until the extrapolation condition is met - double time = starttime; - double tstart = time; - // cosmics are always (?) downwards, we can simplify this logic using that TODO - while(fabs(time-tstart) < TCRV.maxDt() && TCRV.needsExtrapolation(ftraj,tdir) ){ - TimeRange range = tdir == TimeDir::forwards ? TimeRange(time,time+TCRV.step()) : TimeRange(time-TCRV.step(),time); - ktrk.extendTraj(range); - time = tdir == TimeDir::forwards ? range.end() : range.begin(); - } - hadintersection = false; - if (TCRV.intersection().good()){ - hadintersection = true; - // we have a good intersection. Use this to create a Shell material Xing - auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); - // TODO add DS and shielding material - KKCRVXINGPTR crvxingptr = std::make_shared(TCRV.module(), TCRVSID, *kkmat_h->STMaterial(),TCRV.intersection(),reftrajptr,tcrvthick_,TCRV.interTolerance()); - ktrk.addTCRVXing(crvxingptr,tdir); - } - } while(hadintersection); - } - } #endif diff --git a/Mu2eKinKal/inc/KKFit.hh b/Mu2eKinKal/inc/KKFit.hh index 36f5d2bb87..32e47d6aba 100644 --- a/Mu2eKinKal/inc/KKFit.hh +++ b/Mu2eKinKal/inc/KKFit.hh @@ -7,7 +7,7 @@ #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" -#include "Offline/KinKalGeom/inc/KKMaterial.hh" +#include "Offline/KinKalGeom/inc/KKStrawMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKCaloHit.hh" #include "Offline/Mu2eKinKal/inc/KKFitUtilities.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" @@ -96,18 +96,18 @@ namespace mu2e { explicit KKFit(KKFitConfig const& fitconfig); // helper functions used to create components of the fit // Make KKStrawHits from ComboHits - bool makeStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf, + bool makeStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf, KKStrawMaterial const& smat, PTRAJ const& ptraj, ComboHitCollection const& chcol, StrawHitIndexCollection const& strawHitIdxs, KKSTRAWHITCOL& hits, KKSTRAWXINGCOL& exings) const; // regrow KKTrack components from a KalSeed bool regrowComponents(KalSeed const& kseed, ComboHitCollection const& chcol, mu2e::IndexMap const& strawindexmap, - Tracker const& tracker,Calorimeter const& calo, StrawResponse const& strawresponse,BFieldMap const& kkbf, + Tracker const& tracker,Calorimeter const& calo, StrawResponse const& strawresponse,BFieldMap const& kkbf, KKStrawMaterial const& smat, PTRAJPTR& ptraj, KKSTRAWHITCOL& strawhits, KKCALOHITCOL& calohits, KKSTRAWXINGCOL& exings, DOMAINCOL& domains) const; std::shared_ptr caloAxis(CaloCluster const& cluster, Calorimeter const& calo) const; // should come from CaloCluster TODO bool makeCaloHit(CCPtr const& cluster, Calorimeter const& calo, PTRAJ const& pktraj, KKCALOHITCOL& hits) const; // extend a track with a new configuration, optionally searching for and adding hits and straw material void extendTrack(Config const& config, BFieldMap const& kkbf, Tracker const& tracker, - StrawResponse const& strawresponse, ComboHitCollection const& chcol, + StrawResponse const& strawresponse, KKStrawMaterial const& smat, ComboHitCollection const& chcol, Calorimeter const& calo, CCHandle const& cchandle, KKTRK& kktrk) const; // sample the fit at the specificed surfaces @@ -122,9 +122,9 @@ namespace mu2e { auto const& strawHitClusterer() const { return shclusterer_; } private: void fillTrackerInfo(Tracker const& tracker) const; - void addStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf, + void addStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf, KKStrawMaterial const& smat, KKTRK const& kktrk, ComboHitCollection const& chcol, KKSTRAWHITCOL& hits,KKSTRAWXINGCOL& addexings) const; - void addStraws(Tracker const& tracker, KKTRK const& kktrk, KKSTRAWXINGCOL& addexings) const; + void addStraws(Tracker const& tracker, KKStrawMaterial const& smat, KKTRK const& kktrk, KKSTRAWXINGCOL& addexings) const; void addCaloHit(Calorimeter const& calo, KKTRK& kktrk, CCHandle cchandle, KKCALOHITCOL& hits) const; void sampleFit(KKTRK const& kktrk,KalIntersectionCollection& inters) const; // sample fit at the surfaces specified in the config int print_; @@ -212,12 +212,10 @@ namespace mu2e { } } - template bool KKFit::makeStrawHits(Tracker const& tracker,StrawResponse const& strawresponse,BFieldMap const& kkbf, + template bool KKFit::makeStrawHits(Tracker const& tracker,StrawResponse const& strawresponse,BFieldMap const& kkbf, KKStrawMaterial const& smat, PTRAJ const& ptraj, ComboHitCollection const& chcol, StrawHitIndexCollection const& strawHitIdxs, KKSTRAWHITCOL& hits, KKSTRAWXINGCOL& exings) const { unsigned ngood(0); - GeomHandle kkmat_h; - auto const& smat = kkmat_h->strawMaterial(); // loop over the individual straw combo hits for(auto strawidx : strawHitIdxs) { const ComboHit& combohit(chcol.at(strawidx)); @@ -254,10 +252,8 @@ namespace mu2e { template bool KKFit::regrowComponents(KalSeed const& kseed, // primary event input ComboHitCollection const& chcol, mu2e::IndexMap const& strawindexmap, // ancillary event input Tracker const& tracker,Calorimeter const& calo, // geometries - StrawResponse const& strawresponse,BFieldMap const& kkbf, // other conditions + StrawResponse const& strawresponse,BFieldMap const& kkbf, KKStrawMaterial const& smat, // other conditions PTRAJPTR& ptraj, KKSTRAWHITCOL& strawhits, KKCALOHITCOL& calohits, KKSTRAWXINGCOL& exings, DOMAINCOL& domains) const { // return values - GeomHandle kkmat_h; - auto const& smat = kkmat_h->strawMaterial(); unsigned ngood(0), nactive(0), nsactive(0); // loop over the TrkStrawHitSeeds in this KalSeed for(auto const& tshs : kseed.hits()) { @@ -353,14 +349,14 @@ namespace mu2e { } template void KKFit::extendTrack(Config const& exconfig, BFieldMap const& kkbf, Tracker const& tracker, - StrawResponse const& strawresponse, ComboHitCollection const& chcol, + StrawResponse const& strawresponse, KKStrawMaterial const& smat, ComboHitCollection const& chcol, Calorimeter const& calo, CCHandle const& cchandle, KKTRK& kktrk) const { KKSTRAWHITCOL addstrawhits; KKCALOHITCOL addcalohits; KKSTRAWXINGCOL addstrawxings; - if(addhits_)addStrawHits(tracker, strawresponse, kkbf, kktrk, chcol, addstrawhits, addstrawxings ); - if(matcorr_ && addmat_)addStraws(tracker, kktrk, addstrawxings); + if(addhits_)addStrawHits(tracker, strawresponse, kkbf, smat, kktrk, chcol, addstrawhits, addstrawxings ); + if(matcorr_ && addmat_)addStraws(tracker, smat, kktrk, addstrawxings); if(addhits_ && usecalo_ && kktrk.caloHits().size()==0)addCaloHit(calo, kktrk, cchandle, addcalohits); if(print_ > 1){ std::cout << "KKTrk extension adding " @@ -371,10 +367,8 @@ namespace mu2e { kktrk.extendTrack(exconfig,addstrawhits,addstrawxings,addcalohits); } - template void KKFit::addStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf, + template void KKFit::addStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf, KKStrawMaterial const& smat, KKTRK const& kktrk, ComboHitCollection const& chcol, KKSTRAWHITCOL& addhits, KKSTRAWXINGCOL& addexings) const { - GeomHandle kkmat_h; - auto const& smat = kkmat_h->strawMaterial(); auto const& ptraj = kktrk.fitTraj(); // add the buffer to the time range; this defines the search range for new hits TimeRange brange(ptraj.range().begin()-maxStrawHitDt_, ptraj.range().end()+maxStrawHitDt_); @@ -435,11 +429,8 @@ namespace mu2e { } } - template void KKFit::addStraws(Tracker const& tracker, KKTRK const& kktrk, + template void KKFit::addStraws(Tracker const& tracker, KKStrawMaterial const& smat, KKTRK const& kktrk, KKSTRAWXINGCOL& addexings) const { - GeomHandle kkmat_h; - auto const& smat = kkmat_h->strawMaterial(); - // this algorithm assumes the track never hits the same straw twice. That could be violated by reflecting tracks, and could be addressed // by including the time of the Xing as part of its identity. That would slow things down so it remains to be proven it's a problem TODO // build the set of existing straws diff --git a/Mu2eKinKal/inc/KKFitSettings.hh b/Mu2eKinKal/inc/KKFitSettings.hh index d1dd609765..3a881c64f6 100644 --- a/Mu2eKinKal/inc/KKFitSettings.hh +++ b/Mu2eKinKal/inc/KKFitSettings.hh @@ -112,12 +112,10 @@ namespace mu2e { fhicl::Atom btol { Name("BCorrTolerance"), Comment("Tolerance on BField correction momentum fractional accuracy (dimensionless)") }; fhicl::Atom interTol { Name("IntersectionTolerance"), Comment("Tolerance for surface intersections (mm)") }; fhicl::Atom MaxDt { Name("MaxDt"), Comment("Maximum time to extrapolate a fit (ns)") }; - fhicl::Atom MinV { Name("MinV"), Comment("Minimum Y vel to extrapolate a fit") }; fhicl::Atom BackToTracker { Name("BackToTracker"), Comment("Extrapolate reflecting tracks back to the tracker") }; fhicl::Atom ToTrackerEnds { Name("ToTrackerEnds"), Comment("Extrapolate tracks to the tracker ends") }; fhicl::Atom Upstream { Name("Upstream"), Comment("Extrapolate tracks upstream") }; fhicl::Atom ToOPA { Name("ToOPA"), Comment("Test tracks for intersection with the OPA") }; - fhicl::Atom ToCRV { Name("ToCRV"), Comment("Extrapolate tracks to the CRV") }; }; } } diff --git a/Mu2eKinKal/src/CentralHelixFit_module.cc b/Mu2eKinKal/src/CentralHelixFit_module.cc index 90ba0e1ef2..f6cef066ed 100644 --- a/Mu2eKinKal/src/CentralHelixFit_module.cc +++ b/Mu2eKinKal/src/CentralHelixFit_module.cc @@ -52,6 +52,7 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" @@ -241,6 +242,7 @@ namespace mu2e { void CentralHelixFit::produce(art::Event& event ) { GeomHandle calo_h; GeomHandle nominalTracker_h; + GeomHandle kkmat_h; // find current proditions auto const& strawresponse = strawResponse_h_.getPtr(event.id()); auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); @@ -314,7 +316,7 @@ namespace mu2e { strawhits.reserve(strawHitIdxs.size()); KKSTRAWXINGCOL strawxings; strawxings.reserve(strawHitIdxs.size()); - kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, pseedtraj, chcol, strawHitIdxs, strawhits, strawxings); + kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_h->strawMaterial(), pseedtraj, chcol, strawHitIdxs, strawhits, strawxings); // optionally (and if present) add the CaloCluster as a constraint // verify the cluster looks physically reasonable before adding it TODO! Or, let the KKCaloHit updater do it TODO KKCALOHITCOL calohits; @@ -330,7 +332,7 @@ namespace mu2e { // if we have an extension schedule, extend. if(goodfit && exconfig_.schedule().size() > 0) { // std::cout << "EXTENDING TRACK " << event.id() << " " << index << std::endl; - kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, chcol, *calo_h, cc_H, *kktrk ); + kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H, *kktrk ); goodfit = goodFit(*kktrk); } diff --git a/Mu2eKinKal/src/KinematicLineFit_module.cc b/Mu2eKinKal/src/KinematicLineFit_module.cc index acf1594f6c..9f309c0ba8 100644 --- a/Mu2eKinKal/src/KinematicLineFit_module.cc +++ b/Mu2eKinKal/src/KinematicLineFit_module.cc @@ -52,7 +52,7 @@ #include "Offline/Mu2eKinKal/inc/KKBField.hh" #include "Offline/Mu2eKinKal/inc/KKFitUtilities.hh" #include "Offline/Mu2eKinKal/inc/KKShellXing.hh" -#include "Offline/Mu2eKinKal/inc/KKExtrap.hh" +#include "Offline/Mu2eKinKal/inc/ExtrapolateTCRV.hh" // root #include "TH1F.h" #include "TTree.h" @@ -119,12 +119,20 @@ namespace mu2e { fhicl::Atom sampleTBuff { Name("SampleTimeBuffer"), Comment("Time buffer for sample intersections (nsec)") }; }; + // Extrapolation configuration + struct KKExtrapConfig { + fhicl::Atom MaxDt { Name("MaxDt"), Comment("Maximum time to extrapolate a fit") }; + fhicl::Atom MinV { Name("MinV"), Comment("Minimum Y vel to extrapolate a fit") }; + fhicl::Atom ToCRV { Name("ToCRV"), Comment("Extrapolate tracks to the CRV") }; + fhicl::Atom Debug { Name("Debug"), Comment("Debug level"), 0 }; + }; + struct GlobalConfig { fhicl::Table modSettings { Name("ModuleSettings") }; fhicl::Table mu2eSettings { Name("KKFitSettings") }; fhicl::Table fitSettings { Name("FitSettings") }; fhicl::Table extSettings { Name("ExtensionSettings") }; - fhicl::OptionalTable Extrapolation { Name("ExtrapolationSettings") }; + fhicl::OptionalTable Extrapolation { Name("Extrapolation") }; fhicl::OptionalAtom fitDirection { Name("FitDirection"), Comment("Particle direction to fit, either \"upstream\" or \"downstream\"")}; }; @@ -161,9 +169,12 @@ namespace mu2e { double intertol_; // surface intersection tolerance (mm) double sampletbuff_; // simple time buffer; replace this with extrapolation TODO bool sampleinrange_, sampleinbounds_; // require samples to be in range or on surface - std::unique_ptr extrap_; // extrapolation helper + bool extrapolate_, toCRV_; + double maxdt_ = 0.0, btol_ = 0.0, minv_ = 0.0; SurfaceIdCollection ssids_; KinKalGeom::SurfacePairCollection surfacess_to_sample_; // surfaces to sample the fit + int extrapdebug_ = 0; + double tcrvthick_ = 150.0; // CRV sector thickness: should come from geometry service TODO Config config_; // initial fit configuration object Config exconfig_; // extension configuration object }; @@ -181,6 +192,7 @@ namespace mu2e { sampletbuff_(settings().modSettings().sampleTBuff()), sampleinrange_(settings().modSettings().sampleInRange()), sampleinbounds_(settings().modSettings().sampleInBounds()), + extrapolate_(false), toCRV_(false), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), exconfig_(Mu2eKinKal::makeConfig(settings().extSettings())) { @@ -211,9 +223,21 @@ namespace mu2e { for(auto const& sidname : settings().modSettings().sampleSurfaces()) { ssids_.push_back(SurfaceId(sidname,-1)); // match all elements } - // setup extrapolation - if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings()); + // configure extrapolation + if(settings().Extrapolation()){ + extrapolate_ = true; + toCRV_ = settings().Extrapolation()->ToCRV(); + // global configs + maxdt_ = settings().Extrapolation()->MaxDt(); + btol_ = settings().extSettings().btol(); // use the same BField cor. tolerance as in fit extension + minv_ = settings().Extrapolation()->MinV(); + extrapdebug_ = settings().Extrapolation()->Debug(); + } + + if(print_ > 0) std::cout << config_; + + } KinematicLineFit::~KinematicLineFit(){} @@ -235,6 +259,7 @@ namespace mu2e { void KinematicLineFit::produce(art::Event& event ) { GeomHandle calo_h; GeomHandle nominalTracker_h; + GeomHandle kkmat_h; // find current proditions auto const& strawresponse = strawResponse_h_.getPtr(event.id()); @@ -275,7 +300,7 @@ namespace mu2e { KKSTRAWXINGCOL strawxings; strawhits.reserve(strawHitIdxs.size()); strawxings.reserve(strawHitIdxs.size()); - kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, pseedtraj, *chcolptr, strawHitIdxs, strawhits, strawxings); + kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_h->strawMaterial(), pseedtraj, *chcolptr, strawHitIdxs, strawhits, strawxings); //here KKCALOHITCOL calohits; @@ -296,11 +321,11 @@ namespace mu2e { auto kktrk = make_unique(config_,*kkbf_,seedtraj,fpart_,kkfit_.strawHitClusterer(),strawhits,strawxings,calohits,paramconstraints_); auto goodfit = goodFit(*kktrk); if(goodfit && exconfig_.schedule().size() > 0){ - kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, chcol, *calo_h, cc_H, *kktrk ); + kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H, *kktrk ); } goodfit = goodFit(*kktrk); // extrapolate as required - if(goodfit && extrap_) extrap_->extrapolate(*kktrk); + if(goodfit && extrapolate_) extrapolate(*kktrk); bool save = goodFit(*kktrk); if(save || saveall_){ TrkFitFlag fitflag(hptr->status()); @@ -380,5 +405,38 @@ namespace mu2e { } } } + + void KinematicLineFit::extrapolate(KKTRK& ktrk) const { + GeomHandle kkg_h; + GeomHandle kkmat_h; + // extrapolate to the extracted CRV. This function should be migrated to KKExtrap TODO + auto TCRV = ExtrapolateTCRV(maxdt_,btol_,intertol_,minv_,*kkg_h->TCRV(),extrapdebug_); + + auto const& ftraj = ktrk.fitTraj(); + static const SurfaceId TCRVSID("TCRV"); + auto dir0 = ftraj.direction(ftraj.t0()); + TimeDir tdir = (dir0.Y() > 0) ? TimeDir::forwards : TimeDir::backwards; + double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); + bool hadintersection = false; + do { + // iterate until the extrapolation condition is met + double time = starttime; + double tstart = time; + while(fabs(time-tstart) < TCRV.maxDt() && TCRV.needsExtrapolation(ftraj,tdir) ){ + TimeRange range = tdir == TimeDir::forwards ? TimeRange(time,time+TCRV.step()) : TimeRange(time-TCRV.step(),time); + ktrk.extendTraj(range); + time = tdir == TimeDir::forwards ? range.end() : range.begin(); + } + hadintersection = false; + if (TCRV.intersection().good()){ + hadintersection = true; + // we have a good intersection. Use this to create a Shell material Xing + auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); + // TODO add DS and shielding material + KKCRVXINGPTR crvxingptr = std::make_shared(TCRV.module(), TCRVSID, *kkmat_h->STMaterial(),TCRV.intersection(),reftrajptr,tcrvthick_,TCRV.interTolerance()); + ktrk.addTCRVXing(crvxingptr,tdir); + } + } while(hadintersection); + } } DEFINE_ART_MODULE(mu2e::KinematicLineFit) diff --git a/Mu2eKinKal/src/LoopHelixFit_module.cc b/Mu2eKinKal/src/LoopHelixFit_module.cc index 1ecec5900f..fca260082a 100644 --- a/Mu2eKinKal/src/LoopHelixFit_module.cc +++ b/Mu2eKinKal/src/LoopHelixFit_module.cc @@ -52,6 +52,7 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" @@ -285,6 +286,7 @@ namespace mu2e { throw cet::exception("RECO") << "mu2e::LoopHelixFit: Unknown helix propagation direction " << fdir.name(); // geom GeomHandle calo_h; + GeomHandle kkmat_h; // find current proditions auto const& strawresponse = strawResponse_h_.getPtr(event.id()); auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); @@ -329,7 +331,7 @@ namespace mu2e { strawhits.reserve(strawHitIdxs.size()); KKSTRAWXINGCOL strawxings; strawxings.reserve(strawHitIdxs.size()); - if(!kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, pseedtraj, chcol, strawHitIdxs, strawhits, strawxings)) { + if(!kkfit_.makeStrawHits(*tracker, *strawresponse, *kkbf_, kkmat_h->strawMaterial(), pseedtraj, chcol, strawHitIdxs, strawhits, strawxings)) { if(print_>0) printf("[LoopHelixFit::%s] Failed to create a track\n", __func__); return nullptr; } @@ -351,7 +353,7 @@ namespace mu2e { __func__, goodfit, ktrk->fitStatus().chisq_.probability(), ktrk->strawHits().size(), ktrk->caloHits().size()); // if we have an extension schedule, extend. if(goodfit && exconfig_.schedule().size() > 0) { - kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, chcol, *calo_h, cc_H, *ktrk ); + kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H, *ktrk ); goodfit = goodFit(*ktrk,seedtraj); // if finalizing, apply that now. if(goodfit && fconfig_.schedule().size() > 0){ diff --git a/Mu2eKinKal/src/RegrowKinematicLine_module.cc b/Mu2eKinKal/src/RegrowKinematicLine_module.cc index ff02d69854..179331d13c 100644 --- a/Mu2eKinKal/src/RegrowKinematicLine_module.cc +++ b/Mu2eKinKal/src/RegrowKinematicLine_module.cc @@ -54,6 +54,7 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" @@ -178,6 +179,7 @@ namespace mu2e { auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); GeomHandle nominalTracker_h; GeomHandle calo_h; + GeomHandle kkmat_h; // find input event data auto kseed_H = event.getValidHandle(kseedcol_T_); @@ -212,7 +214,7 @@ namespace mu2e { DOMAINCOL domains; // create the trajectory. This is done here to be strongly typed auto goodhits = kkfit_.regrowComponents(kseed, chcol, indexmap, - *tracker,*calo_h,*strawresponse,*kkbf_, + *tracker,*calo_h,*strawresponse,*kkbf_, kkmat_h->strawMaterial(), trajptr, strawhits, calohits, strawxings, domains); if(debug_ > 0){ std::cout << "Regrew " << strawhits.size() << " straw hits, " << strawxings.size() << " straw xings, " << calohits.size() << " CaloHits and " << domains.size() << " domains, status = " << goodhits << std::endl; diff --git a/Mu2eKinKal/src/RegrowLoopHelix_module.cc b/Mu2eKinKal/src/RegrowLoopHelix_module.cc index 085ed2166d..a4584dfe52 100644 --- a/Mu2eKinKal/src/RegrowLoopHelix_module.cc +++ b/Mu2eKinKal/src/RegrowLoopHelix_module.cc @@ -54,6 +54,7 @@ #include "Offline/Mu2eKinKal/inc/KKFit.hh" #include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" +#include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" #include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" @@ -78,6 +79,7 @@ namespace mu2e { using Mu2eKinKal::KKFinalConfig; using KKFitConfig = Mu2eKinKal::KKFitConfig; using KKModuleConfig = Mu2eKinKal::KKModuleConfig; + using KKMaterialConfig = KKMaterial::Config; using SDIS = std::set; using Name = fhicl::Name; @@ -130,6 +132,8 @@ namespace mu2e { using DOMAINPTR = std::shared_ptr; using DOMAINCOL = std::set; + using KKMaterialConfig = KKMaterial::Config; + explicit RegrowLoopHelix(const Parameters& settings); void beginRun(art::Run& run) override; void produce(art::Event& event) override; @@ -185,6 +189,7 @@ namespace mu2e { auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); GeomHandle nominalTracker_h; GeomHandle calo_h; + GeomHandle kkmat_h; // find input event data auto kseed_H = event.getValidHandle(kseedcol_T_); const auto& kseedcol = *kseed_H; @@ -219,7 +224,7 @@ namespace mu2e { DOMAINCOL domains; // create the trajectory. This is done here to be strongly typed auto goodhits = kkfit_.regrowComponents(kseed, chcol, indexmap, - *tracker,*calo_h,*strawresponse,*kkbf_, + *tracker,*calo_h,*strawresponse,*kkbf_, kkmat_h->strawMaterial(), trajptr, strawhits, calohits, strawxings, domains); if(debug_ > 1){ std::cout << "Regrew " << strawhits.size() << " straw hits, " << strawxings.size() << " straw xings, " << calohits.size() << " CaloHits and " << domains.size() << " domains, status = " << goodhits << std::endl; @@ -238,7 +243,7 @@ namespace mu2e { auto ktrk = std::make_unique(config_,*kkbf_,kseed.particle(),trajptr,strawhits,strawxings,calohits,domains); if(ktrk && ktrk->fitStatus().usable()){ if(debug_ > 0) std::cout << "RegrowLoopHelix: successful track refit" << std::endl; - if(extend_)kkfit_.extendTrack(config_,*kkbf_, *tracker,*strawresponse, chcol, *calo_h, cc_H , *ktrk ); + if(extend_)kkfit_.extendTrack(config_,*kkbf_, *tracker,*strawresponse, kkmat_h->strawMaterial(), chcol, *calo_h, cc_H , *ktrk ); if(ktrk->fitStatus().usable()){ // extrapolate as requested if(extrap_)extrap_->extrapolate(*ktrk);