diff --git a/.muse b/.muse index 07282d58ce..dc9157ae41 100644 --- a/.muse +++ b/.muse @@ -1,5 +1,5 @@ # prefer to build with this environment -ENVSET p084 +ENVSET p086 # add Offline/bin to path PATH bin # recent commits can take enforce these flags diff --git a/CMakeLists.txt b/CMakeLists.txt index c29badb003..19beee1b27 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -30,6 +30,7 @@ cet_set_compiler_flags(DIAGS VIGILANT -Wno-unused-parameter -Wno-non-virtual-dtor -Wno-extra + -Wreorder ) if(DEFINED Offline_UPS_QUALIFIER_STRING) @@ -54,7 +55,7 @@ endif() message("Creating link from ${CMAKE_CURRENT_SOURCE_DIR} to ${CMAKE_CURRENT_SOURCE_DIR}/Offline to satisfy includes") file(CREATE_LINK ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/Offline SYMBOLIC) - + find_package(art_root_io REQUIRED EXPORT) find_package(GSL REQUIRED EXPORT) find_package(PostgreSQL REQUIRED EXPORT) @@ -64,8 +65,8 @@ find_package(XercesC REQUIRED EXPORT) find_package(BLAS REQUIRED EXPORT) find_package(artdaq-core-mu2e REQUIRED EXPORT) find_package(TRACE REQUIRED EXPORT) -if( ${WITH_G4} ) - message("--> ADDING G4 LIBS") +if( ${WITH_G4} ) + message("--> ADDING G4 LIBS") find_package(Geant4 REQUIRED EXPORT) # TODO: Find or implement FindCRY.cmake @@ -142,7 +143,7 @@ add_subdirectory(MCDataProducts) add_subdirectory(MECOStyleProtonAbsorberGeom) add_subdirectory(Mu2e) add_subdirectory(Mu2eBTrk) -if( ${WITH_G4} ) +if( ${WITH_G4} ) message("---> ADDING G4 DIRS") add_subdirectory(EventGenerator) add_subdirectory(Mu2eG4) diff --git a/DataProducts/inc/IndexMap.hh b/DataProducts/inc/IndexMap.hh index 673ca76c0e..dc094c7900 100644 --- a/DataProducts/inc/IndexMap.hh +++ b/DataProducts/inc/IndexMap.hh @@ -48,6 +48,7 @@ namespace mu2e { ost << std::endl; } } + auto const& map() const { return _theMap; } private: std::map _theMap; // FullIndex is first, CondensedIndex is second diff --git a/EventDisplay/src/DataInterface.cc b/EventDisplay/src/DataInterface.cc index a1cb6f7b7c..7fed8e16e5 100644 --- a/EventDisplay/src/DataInterface.cc +++ b/EventDisplay/src/DataInterface.cc @@ -1003,7 +1003,7 @@ void DataInterface::fillEvent(boost::shared_ptr const &contentS { const mu2e::TrkStrawHitSeed &hit = hits.at(j); int sid = hit.strawId().asUint16(); - double time = hit.hitTime(); + double time = hit.time(); std::map >::iterator straw=_straws.find(sid); if(straw!=_straws.end() && !std::isnan(time)) diff --git a/Mu2eKinKal/CMakeLists.txt b/Mu2eKinKal/CMakeLists.txt index 038e0998a1..018e12407b 100644 --- a/Mu2eKinKal/CMakeLists.txt +++ b/Mu2eKinKal/CMakeLists.txt @@ -72,7 +72,19 @@ cet_build_plugin(LoopHelixFit art::module Offline::TrkReco ) +cet_build_plugin(RegrowLoopHelix art::module + REG_SOURCE src/RegrowLoopHelix_module.cc + LIBRARIES REG + Offline::Mu2eKinKal + Offline::TrkReco +) +cet_build_plugin(RegrowKinematicLine art::module + REG_SOURCE src/RegrowKinematicLine_module.cc + LIBRARIES REG + Offline::Mu2eKinKal + Offline::TrkReco +) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/data/TrainBkgFinal.dat ${CURRENT_BINARY_DIR} data/TrainBkgFinal.dat COPYONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/data/TrainBkgSeed.dat ${CURRENT_BINARY_DIR} data/TrainBkgSeed.dat COPYONLY) diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 09086e8cf4..592eb822f7 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -43,7 +43,12 @@ Mu2eKinKal : { MaxStrawDOCA : 5.0 # mm MaxStrawDOCAConsistency : 1.0 # units of chi MaxStrawUposBuffer : 0.0 # units of mm + IntersectionTolerance : 0.1 # tolerance for intersections (mm) + SampleInRange : true # require sample be in time range + SampleInBounds : true # require sample be in surface bounds + SampleSurfaces : [] # specific to the fit type SaveHitCalibInfo : false + SaveDomains : false } CHSEEDFIT: { @@ -185,7 +190,7 @@ Mu2eKinKal : { ] } - SEEDFIT: { + LHSEEDFIT: { PrintLevel : 0 MinNDOF : 1 MaxNIter : 10 @@ -194,7 +199,7 @@ Mu2eKinKal : { DivergenceDeltaChisq : 10.0 DivergenceDeltaParams : 1e6 DivergenceGap : 10 # mm - BFieldCorrection : false + BFieldCorrection : true BCorrTolerance : 1e-2 # momemntum fraction ProcessEnds : false MetaIterationSettings : [ @@ -234,7 +239,7 @@ Mu2eKinKal : { StrawXingUpdaterSettings : [] } - SEEDEXT: { + LHSEEDEXT: { PrintLevel : 0 MinNDOF : 1 MaxNIter : 10 @@ -346,24 +351,45 @@ Mu2eKinKal : { CaloClusterCollection : "CaloClusterMaker" } - LOOPHELIX : { + # Minimal common configuration for regrowing KalSeeds + REGROW : { + PrintLevel : 0 + MinNDOF : 1 + MaxNIter : 10 + Deweight : 1.0e6 + ConvergenceDeltaChisq : 0.1 + DivergenceDeltaChisq : 10.0 + DivergenceDeltaParams : 1e6 + DivergenceGap : 100 # mm + ProcessEnds : true + BCorrTolerance : 1e-4 + MetaIterationSettings : [ [ 0.0, "" ] ] + CADSHUSettings : [ ] + BkgANNSHUSettings : [ ] + DriftANNSHUSettings : [ ] + Chi2SHUSettings : [ ] + StrawXingUpdaterSettings : [ [-1.0, 1.0, false, 0 ] ] + } + + LHHELIX : { SeedErrors : [5.0, 5.0, 5.0, 5.0, 0.02, 5.0] # R(mm), Lambda(mm), Cx(mm), Cy(mm), phi0, t0 (ns) SeedFlags : [ "HelixOK" ] - IntersectionTolerance : 0.1 # tolerance for intersections (mm) - SampleInRange : true # require sample be in time range - SampleInBounds : true # require sample be in surface bounds } - LOOPEXTRAPOLATION : { + LHDRIFTXTRAP : { + IntersectionTolerance : 0.1 MaxDt : 200.0 # (ns) + BCorrTolerance : 1e-4 # momemntum fraction ToTrackerEnds : true Upstream : true BackToTracker : false ToOPA : true } - SEEDEXTRAPOLATION : { + LHSEEDXTRAP : { + IntersectionTolerance : 0.1 MaxDt : 200.0 # (ns) + BCorrTolerance : 1e-2 # momemntum fraction ToTrackerEnds : true Upstream : false BackToTracker : false @@ -416,16 +442,16 @@ Mu2eKinKal : { MaterialSettings : @local::Mu2eKinKal.MAT KKFitSettings: { @table::Mu2eKinKal.KKFIT + SaveDomains : false SaveTrajectory : Detector } - FitSettings : @local::Mu2eKinKal.SEEDFIT - ExtensionSettings : @local::Mu2eKinKal.SEEDEXT + FitSettings : @local::Mu2eKinKal.LHSEEDFIT + ExtensionSettings : @local::Mu2eKinKal.LHSEEDEXT ModuleSettings : { - @table::Mu2eKinKal.LOOPHELIX + @table::Mu2eKinKal.LHHELIX @table::Mu2eKinKal.KKPrecursors - SampleSurfaces : [] } - Extrapolation : @local::Mu2eKinKal.SEEDEXTRAPOLATION + ExtrapolationSettings : @local::Mu2eKinKal.LHSEEDXTRAP UsePDGCharge: false HelixMask: { MinHelixMom : 0 @@ -437,17 +463,20 @@ Mu2eKinKal : { 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 # save trajectories in the Detector region + SaveDomains : true SaveTrajectory : Full } - FitSettings : @local::Mu2eKinKal.SEEDFIT + FitSettings : { + @table::Mu2eKinKal.LHSEEDFIT + } ExtensionSettings : @local::Mu2eKinKal.LHDRIFTEXT ModuleSettings : { - @table::Mu2eKinKal.LOOPHELIX + @table::Mu2eKinKal.LHHELIX @table::Mu2eKinKal.KKPrecursors - SampleSurfaces : ["ST_Outer","ST_Front","ST_Back"] # these are additional surfaces; surfaces used in extrapolation are also sampled } - Extrapolation : @local::Mu2eKinKal.LOOPEXTRAPOLATION + ExtrapolationSettings : @local::Mu2eKinKal.LHDRIFTXTRAP UsePDGCharge: false HelixMask: { MinHelixMom : 0 @@ -459,7 +488,10 @@ Mu2eKinKal : { MaterialSettings : @local::Mu2eKinKal.MAT KKFitSettings: @local::Mu2eKinKal.KKFIT FitSettings : @local::Mu2eKinKal.CHSEEDFIT - ExtensionSettings : @local::Mu2eKinKal.CHSEEDEXT + ExtensionSettings : { + @table::Mu2eKinKal.CHSEEDEXT + BFieldCorrection : false + } ModuleSettings : { @table::Mu2eKinKal.KINEMATICLINE @table::Mu2eKinKal.KKPrecursors @@ -478,6 +510,7 @@ Mu2eKinKal : { MaxCaloClusterDt: 8 # save the full trajectory SaveTrajectory: Full + SaveDomains : false # meaningless for Linefit } FitSettings : @local::Mu2eKinKal.CHSEEDFIT ExtensionSettings : @local::Mu2eKinKal.CHDRIFTEXT @@ -494,6 +527,7 @@ Mu2eKinKal : { KKFitSettings: { @table::Mu2eKinKal.KKFIT SaveTrajectory : T0 + SaveDomains : false } FitSettings : @local::Mu2eKinKal.CHSEEDFIT ExtensionSettings : @local::Mu2eKinKal.CHSEEDEXT @@ -509,6 +543,7 @@ Mu2eKinKal : { KKFitSettings: { @table::Mu2eKinKal.KKFIT SaveTrajectory : Full + SaveDomains : true } FitSettings : @local::Mu2eKinKal.CHSEEDFIT ExtensionSettings : @local::Mu2eKinKal.CHDRIFTEXT @@ -538,15 +573,11 @@ Mu2eKinKal : { KKUmu: @local::Mu2eKinKal.LHDriftFit } } -# Extrapolate upstream fits back to the tracker entrance -Mu2eKinKal.producers.KKUe.Extrapolation.BackToTracker : true -Mu2eKinKal.producers.KKUmu.Extrapolation.BackToTracker : true - +# update directions and particles as needed Mu2eKinKal.producers.KKCHSeedFitmuM.ModuleSettings.FitParticle : @local::Particle.muminus Mu2eKinKal.producers.KKCHSeedFitmuP.ModuleSettings.FitParticle : @local::Particle.muplus Mu2eKinKal.producers.KKCHmu.ModuleSettings.FitParticle : @local::Particle.muminus # charge floats in the fit, this just determines the mass -Mu2eKinKal.producers.KLSeedFit.ExtensionSettings.BFieldCorrection : false Mu2eKinKal.producers.KLSeedFit.ModuleSettings.FitParticle : @local::Particle.muminus Mu2eKinKal.producers.KKLine.ModuleSettings.FitParticle : @local::Particle.muminus @@ -563,17 +594,15 @@ Mu2eKinKal.producers.KKUmuSeedFit.FitDirection : @local::FitDir.downstream # save trajectories in the Detector region for downstream fits, just the T0 segment for the rest Mu2eKinKal.producers.KKDe.ModuleSettings.FitParticle : @local::Particle.eminus Mu2eKinKal.producers.KKDe.FitDirection : @local::FitDir.downstream -Mu2eKinKal.producers.KKDe.KKFitSettings.SaveTrajectory : "Detector" Mu2eKinKal.producers.KKDmu.ModuleSettings.FitParticle : @local::Particle.muminus Mu2eKinKal.producers.KKDmu.FitDirection : @local::FitDir.downstream -Mu2eKinKal.producers.KKDmu.KKFitSettings.SaveTrajectory : "Detector" Mu2eKinKal.producers.KKUe.ModuleSettings.FitParticle : @local::Particle.eminus Mu2eKinKal.producers.KKUe.FitDirection : @local::FitDir.upstream -Mu2eKinKal.producers.KKUe.KKFitSettings.SaveTrajectory : "Detector" Mu2eKinKal.producers.KKUmu.ModuleSettings.FitParticle : @local::Particle.muminus Mu2eKinKal.producers.KKUmu.FitDirection : @local::FitDir.upstream -Mu2eKinKal.producers.KKUmu.KKFitSettings.SaveTrajectory : "Detector" -# extrapolate upstream fits back to the tracker -physics.producers.KKUmu.Extrapolation.BackToTracker : true -physics.producers.KKUe.Extrapolation.BackToTracker : true +# propagate upstream fits back to the tracker +Mu2eKinKal.producers.KKUe.ExtrapolationSettings.BackToTracker : true +Mu2eKinKal.producers.KKUmu.ExtrapolationSettings.BackToTracker : true + +# END_PROLOG diff --git a/Mu2eKinKal/inc/Chi2SHU_updateCluster.hh b/Mu2eKinKal/inc/Chi2SHU_updateCluster.hh index 32d0f1c3b1..a05cd4b5fc 100644 --- a/Mu2eKinKal/inc/Chi2SHU_updateCluster.hh +++ b/Mu2eKinKal/inc/Chi2SHU_updateCluster.hh @@ -72,7 +72,7 @@ namespace mu2e { pvar = resid.parameterVariance(); } // Use the unbiased residual to compute the chisq - Residual uresid(uresidval,resid.variance(),pvar,resid.active(),resid.dRdP()); + Residual uresid(uresidval,resid.variance(),pvar,resid.dRdP(),resid.active()); chisq += uresid.chisq(); ++ndof; } diff --git a/Mu2eKinKal/inc/ExtrapolateIPA.hh b/Mu2eKinKal/inc/ExtrapolateIPA.hh index bd72b26228..b1462569b5 100644 --- a/Mu2eKinKal/inc/ExtrapolateIPA.hh +++ b/Mu2eKinKal/inc/ExtrapolateIPA.hh @@ -57,7 +57,7 @@ namespace mu2e { if(ktraj.range().range() <= epsilon) return true; // keep going if the step is very small auto stime = tdir == TimeDir::forwards ? ktraj.range().begin()+epsilon : ktraj.range().end()-epsilon; auto etime = tdir == TimeDir::forwards ? ktraj.range().end() : ktraj.range().begin(); - auto vel = ktraj.speed(stime)*ktraj.axis(stime).direction();// use helix axis to define the velocity + auto vel = ktraj.speed(stime)*ktraj.linearize(stime).direction();// use hlinear approximation auto spos = ktraj.position3(stime); auto epos = ktraj.position3(etime); double zvel = vel.Z()*timeDirSign(tdir); // sign by extrapolation direction diff --git a/Mu2eKinKal/inc/KKCaloHit.hh b/Mu2eKinKal/inc/KKCaloHit.hh index f80544b5cb..d56185cfa8 100644 --- a/Mu2eKinKal/inc/KKCaloHit.hh +++ b/Mu2eKinKal/inc/KKCaloHit.hh @@ -26,6 +26,29 @@ namespace mu2e { using CA = KinKal::ClosestApproach; using HIT = KinKal::Hit; using KTRAJPTR = std::shared_ptr; + // clone op for reinstantiation + KKCaloHit(KKCaloHit const& rhs): + caloCluster_(rhs.caloCluster()), + axis_(rhs.sensorAxis()), + tvar_(rhs.timeVariance()), + wvar_(rhs.widthVariance()), + ca_( + rhs.closestApproach().particleTraj(), + axis_, + rhs.closestApproach().hint(), + rhs.closestApproach().precision() + ), + rresid_(rhs.timeResidual()){ + /**/ + }; + std::shared_ptr< KinKal::Hit > clone(CloneContext& context) const override{ + auto rv = std::make_shared< KKCaloHit >(*this); + auto ca = rv->closestApproach(); + auto trajectory = std::make_shared(ca.particleTraj()); + ca.setTrajectory(trajectory); + rv->setClosestApproach(ca); + return rv; + }; // Hit interface overrrides unsigned nResid() const override { return 1; } // 1 time residual Residual const& refResidual(unsigned ires=0) const override; @@ -54,6 +77,8 @@ namespace mu2e { double wvar_; // variance in transverse position of the sensor/measurement in mm. Assumes cylindrical error, could be more general CA ca_; // reference time and distance of closest approach to the axis Residual rresid_; // residual WRT most recent reference parameters + // clone support + void setClosestApproach(const CA& ca){ ca_ = ca; } }; template KKCaloHit::KKCaloHit(CCPtr caloCluster, PCA const& pca, double tvar, double wvar) : caloCluster_(caloCluster), axis_(pca.sensorTraj()), tvar_(tvar), wvar_(wvar), @@ -73,7 +98,7 @@ namespace mu2e { // if(ca_.usable()) tphint = CAHint(ca_.particleToca(),ca_.sensorToca()); PCA pca(ptraj,axis_,tphint,precision()); ca_ = pca.localClosestApproach(); - if(!ca_.usable())rresid_ = Residual(rresid_.value(),rresid_.variance(),0.0,false,rresid_.dRdP()); + rresid_ = Residual(rresid_.value(),rresid_.variance(),0.0,rresid_.dRdP(),ca_.usable()); } template void KKCaloHit::updateState(MetaIterConfig const& miconfig,bool first) { @@ -112,9 +137,9 @@ namespace mu2e { double invvar2 = std::max(s2*sint2/(cost2*wvar_), 12/(ldt*ldt)); totvar += 1.0/invvar2; } - rresid_ = Residual(ca_.deltaT(),totvar,0.0,true,ca_.dTdP()); + rresid_ = Residual(ca_.deltaT(),totvar,0.0,ca_.dTdP()); } else { - rresid_ = Residual(rresid_.value(),rresid_.variance(),0.0,false,rresid_.dRdP()); + rresid_ = Residual(rresid_.value(),rresid_.variance(),0.0,rresid_.dRdP(),false); } // finally update the weight this->updateWeight(miconfig); diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh new file mode 100644 index 0000000000..c475b04e1f --- /dev/null +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -0,0 +1,279 @@ +#ifndef Mu2eKinKal_KKExtrap_hh +#define Mu2eKinKal_KKExtrap_hh +// +// Helper class for extrapolating fits through materials and surfaces. +// original author: Dave Brown (LBN), 9/8/2025 +// +#include "Offline/Mu2eKinKal/inc/KKFitSettings.hh" +#include "Offline/DataProducts/inc/SurfaceId.hh" +#include "Offline/KinKalGeom/inc/SurfaceMap.hh" +#include "Offline/KinKalGeom/inc/Tracker.hh" +#include "Offline/Mu2eKinKal/inc/KKTrack.hh" +#include "Offline/Mu2eKinKal/inc/ExtrapolateToZ.hh" +#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 "KinKal/Geometry/ParticleTrajectoryIntersect.hh" +namespace mu2e { + using KinKal::VEC3; + using KinKal::TimeDir; + using SurfacePtr = std::shared_ptr; + using CylPtr = std::shared_ptr; + using DiskPtr = std::shared_ptr; + using FruPtr = std::shared_ptr; + using AnnPtr = std::shared_ptr; + using Mu2eKinKal::KKExtrapConfig; + + class KKExtrap { + public: + explicit KKExtrap(KKExtrapConfig const& exconfig,KKMaterial const& kkmat); + // extrapolation functions; these are templated on the type of trajectory + template void extrapolate(KKTrack& ktrk) const; + template void toTrackerEnds(KKTrack& ktrk) const; + template bool extrapolateIPA(KKTrack& ktrk,TimeDir trkdir) const; + 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 void toOPA(KKTrack& ktrk, double tstart, TimeDir tdir) const; + + private: + int debug_; + double btol_, intertol_, maxdt_; + SurfaceMap smap_; + KKMaterial const& kkmat_; + AnnPtr tsdaptr_; + DiskPtr trkfrontptr_, trkmidptr_, trkbackptr_; + FruPtr opaptr_; + bool backToTracker_, toOPA_, toTrackerEnds_, upstream_; + ExtrapolateToZ TSDA_, trackerFront_, trackerBack_; // extrapolation predicate based on Z values + ExtrapolateIPA extrapIPA_; // extrapolation to intersections with the IPA + ExtrapolateST extrapST_; // extrapolation to intersections with the ST + 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) : + debug_(extrapconfig.Debug()), + btol_(extrapconfig.btol()), + intertol_(extrapconfig.interTol()), + maxdt_(extrapconfig.MaxDt()), + kkmat_(kkmat), + tsdaptr_(smap_.DS().upstreamAbsorberPtr()), + trkfrontptr_(smap_.tracker().frontPtr()), + trkmidptr_(smap_.tracker().middlePtr()), + trkbackptr_(smap_.tracker().backPtr()), + opaptr_(smap_.DS().outerProtonAbsorberPtr()), + backToTracker_(extrapconfig.BackToTracker()), + toOPA_(extrapconfig.ToOPA()), + toTrackerEnds_(extrapconfig.ToTrackerEnds()), + upstream_(extrapconfig.Upstream()), + TSDA_(maxdt_,btol_,smap_.DS().upstreamAbsorber().center().Z(),debug_), + trackerFront_(maxdt_,btol_,smap_.tracker().front().center().Z(),debug_), + trackerBack_(maxdt_,btol_,smap_.tracker().back().center().Z(),debug_), + extrapIPA_(maxdt_,btol_,intertol_,smap_.DS().innerProtonAbsorberPtr(),debug_), + extrapST_(maxdt_,btol_,intertol_,smap_.ST(),debug_) + { + if(debug_ > 0)std::cout << "IPA limits z " << extrapIPA_.zmin() << " " << extrapIPA_.zmax() << std::endl; + if(debug_ > 0)std::cout << "ST limits z " << extrapST_.zmin() << " " << extrapST_.zmax() << " r " << extrapST_.rmin() << " " << extrapST_.rmax() << std::endl; + } + + template void KKExtrap::extrapolate(KKTrack& ktrk) const { + // define the time direction according to the fit direction inside the tracker + auto const& ftraj = ktrk.fitTraj(); + if(toTrackerEnds_)toTrackerEnds(ktrk); + if(upstream_){ + auto dir0 = ftraj.direction(ftraj.t0()); + TimeDir tdir = (dir0.Z() > 0) ? TimeDir::backwards : TimeDir::forwards; + double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); + // extrapolate through the IPA in this direction. + bool exitsIPA = extrapolateIPA(ktrk,tdir); + if(exitsIPA){ // if it exits out the back, extrapolate through the target + bool exitsST = extrapolateST(ktrk,tdir); + if(exitsST) { // if it exits out the back, extrapolate to the TSDA (DS rear absorber) + bool hitTSDA = extrapolateTSDA(ktrk,tdir); + // if we hit the TSDA we are done. Otherwise if we reflected, go back through the ST + if(!hitTSDA){ // reflection upstream of the target: go back through the target + extrapolateST(ktrk,tdir); + if(backToTracker_){ // optionally extrapolate back through the IPA, then to the tracker entrance + extrapolateIPA(ktrk,tdir); + extrapolateTracker(ktrk,tdir); + } + } + } else { // reflection inside the ST; extrapolate back through the IPA, then to the tracker entrance + if(backToTracker_){ + extrapolateIPA(ktrk,tdir); + extrapolateTracker(ktrk,tdir); + } + } + } else { // reflection inside the IPA; extrapolate back through the IPA, then to the tracker entrance + if(backToTracker_)ktrk.extrapolate(tdir,trackerFront_); + } + // optionally test for intersection with the OPA + if(toOPA_)toOPA(ktrk,starttime,tdir); + } + } + + template void KKExtrap::toTrackerEnds(KKTrack& ktrk) const { + // 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 volumei + auto const& ftraj = ktrk.fitTraj(); + auto dir0 = ftraj.direction(ftraj.t0()); + TimeDir fronttdir = (dir0.Z() > 0) ? TimeDir::backwards : TimeDir::forwards; + TimeDir backtdir = (dir0.Z() > 0) ? TimeDir::forwards : TimeDir::backwards; + auto tofront = ktrk.extrapolate(fronttdir,trackerFront_); + auto toback = ktrk.extrapolate(backtdir,trackerBack_); + // record the standard tracker intersections + static const SurfaceId tt_front("TT_Front"); + static const SurfaceId tt_mid("TT_Mid"); + static const SurfaceId tt_back("TT_Back"); + + // start with the middle + auto midinter = KinKal::intersect(ftraj,*trkmidptr_,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,*trkfrontptr_,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,*trkfrontptr_,frange,intertol_,fronttdir); + } + if(frontinter.good()) ktrk.addIntersection(tt_front,frontinter); + } + if(toback){ + // 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,*trkbackptr_,brange,intertol_,backtdir); + if(backinter.good())ktrk.addIntersection(tt_back,backinter); + } + } + + template bool KKExtrap::extrapolateIPA(KKTrack& ktrk,TimeDir tdir) const { + using KKIPAXING = KKShellXing; + using KKIPAXINGPTR = std::shared_ptr; + + if(extrapIPA_.debug() > 0)std::cout << "extrapolating to IPA " << std::endl; + // extraplate the fit through the IPA. This will add material effects for each intersection. It will continue till the + // track exits the IPA + extrapIPA_.reset(); + auto const& ftraj = ktrk.fitTraj(); + static const SurfaceId IPASID("IPA"); + double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); + auto startdir = ftraj.direction(starttime); + do { + ktrk.extrapolate(tdir,extrapIPA_); + 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 = smap_.DS().innerProtonAbsorberPtr(); + KKIPAXINGPTR ipaxingptr = std::make_shared(IPA,IPASID,*kkmat_.IPAMaterial(),extrapIPA_.intersection(),reftrajptr,ipathick_,extrapIPA_.interTolerance()); + if(extrapIPA_.debug() > 0){ + double dmom, paramomvar, perpmomvar; + ipaxingptr->materialEffects(dmom,paramomvar,perpmomvar); + std::cout << "IPA Xing dmom " << dmom << " para momsig " << sqrt(paramomvar) << " perp momsig " << sqrt(perpmomvar) << std::endl; + std::cout << " before append mom = " << reftrajptr->momentum(); + } + ktrk.addIPAXing(ipaxingptr,tdir); + if(extrapIPA_.debug() > 0){ + auto const& newtrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); + std::cout << " after append mom = " << newtrajptr->momentum() << std::endl; + } + } + } while(extrapIPA_.intersection().good()); + // check if the particle exited in the same physical direction or not (reflection) + double endtime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); + auto enddir = ftraj.direction(endtime); + if(enddir.Z() * startdir.Z() > 0.0){ + return true; + } + return false; + } + + template bool KKExtrap::extrapolateST(KKTrack& ktrk,TimeDir tdir) const { + using KKSTXING = KKShellXing; + using KKSTXINGPTR = std::shared_ptr; + + // 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 + extrapST_.reset(); + auto const& ftraj = ktrk.fitTraj(); + double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); + auto startdir = ftraj.direction(starttime); + if(extrapST_.debug() > 0)std::cout << "extrapolating to ST " << std::endl; + do { + ktrk.extrapolate(tdir,extrapST_); + 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()); + if(extrapST_.debug() > 0){ + double dmom, paramomvar, perpmomvar; + stxingptr->materialEffects(dmom,paramomvar,perpmomvar); + std::cout << "ST Xing dmom " << dmom << " para momsig " << sqrt(paramomvar) << " perp momsig " << sqrt(perpmomvar) << std::endl; + std::cout << " before append mom = " << reftrajptr->momentum(); + } + // Add the xing. This truncates the fit + ktrk.addSTXing(stxingptr,tdir); + if(extrapST_.debug() > 0){ + auto const& newtrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); + std::cout << " after append mom = " << newtrajptr->momentum() << std::endl; + } + } + } while(extrapST_.intersection().good()); + // check if the particle exited in the same physical direction or not (reflection) + double endtime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); + auto enddir = ftraj.direction(endtime); + if(enddir.Z() * startdir.Z() > 0.0){ + return true; + } + return false; + } + + template bool KKExtrap::extrapolateTracker(KKTrack& ktrk,TimeDir tdir) const { + 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,*trkfrontptr_,ktraj.range(),intertol_,tdir); + if(trkfrontinter.onsurface_){ // dont worry about bounds here + ktrk.addIntersection(TrackerSID,trkfrontinter); + return true; + } + return false; + } + + template bool KKExtrap::extrapolateTSDA(KKTrack& ktrk,TimeDir tdir) const { + if(TSDA_.debug() > 0)std::cout << "extrapolating to TSDA " << std::endl; + auto const& ftraj = ktrk.fitTraj(); + static const SurfaceId TSDASID("TSDA"); + ktrk.extrapolate(tdir,TSDA_); + // if we reflected we're done. Otherwize, save the TSDA intersection + double tend = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); + auto epos = ftraj.position3(tend); + bool retval = epos.Z() < TSDA_.zVal(); + if(retval){ + auto const& ktraj = tdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); + auto tsdainter = KinKal::intersect(ftraj,*tsdaptr_,ktraj.range(),intertol_,tdir); + if(tsdainter.onsurface_)ktrk.addIntersection(TSDASID,tsdainter); + } + return retval; + } + + template void KKExtrap::toOPA(KKTrack& ktrk, double tstart, TimeDir tdir) const { + 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,*opaptr_,trange,intertol_,tdir); + if(opainter.good()){ + ktrk.addIntersection(OPASID,opainter); + } + } +} +#endif diff --git a/Mu2eKinKal/inc/KKFit.hh b/Mu2eKinKal/inc/KKFit.hh index 2ae93e9e20..68943a7614 100644 --- a/Mu2eKinKal/inc/KKFit.hh +++ b/Mu2eKinKal/inc/KKFit.hh @@ -19,6 +19,7 @@ #include "Offline/TrackerConditions/inc/StrawResponse.hh" #include "Offline/DataProducts/inc/PDGCode.hh" #include "Offline/DataProducts/inc/StrawIdMask.hh" +#include "Offline/DataProducts/inc/IndexMap.hh" #include "Offline/TrackerGeom/inc/Tracker.hh" #include "Offline/CalorimeterGeom/inc/Calorimeter.hh" #include "Offline/RecoDataProducts/inc/KalSeed.hh" @@ -29,11 +30,14 @@ #include "Offline/RecoDataProducts/inc/StrawHitIndex.hh" #include "Offline/RecoDataProducts/inc/KalSeedAssns.hh" #include "Offline/RecoDataProducts/inc/KalIntersection.hh" +#include "Offline/DataProducts/inc/SurfaceId.hh" +#include "Offline/KinKalGeom/inc/SurfaceMap.hh" // geometry #include "Offline/KinKalGeom/inc/Tracker.hh" // KinKal includes #include "KinKal/Fit/Status.hh" #include "KinKal/Fit/Config.hh" +#include "KinKal/Fit/Domain.hh" #include "KinKal/Trajectory/ParticleTrajectory.hh" #include "KinKal/Trajectory/PiecewiseClosestApproach.hh" #include "KinKal/Trajectory/SensorLine.hh" @@ -57,6 +61,7 @@ namespace mu2e { using KKTRK = KKTrack; using KKTRKPTR = std::unique_ptr; using PTRAJ = KinKal::ParticleTrajectory; + using PTRAJPTR = std::unique_ptr; using PCA = KinKal::PiecewiseClosestApproach; using TCA = KinKal::ClosestApproach; using KKHIT = KinKal::Measurement; @@ -81,13 +86,20 @@ namespace mu2e { using EXING = KinKal::ElementXing; using EXINGPTR = std::shared_ptr; using EXINGCOL = std::vector; + using DOMAINPTR = std::shared_ptr; + using DOMAINCOL = std::set; enum SaveTraj {none=0, full, detector, t0seg}; // construct from fit configuration objects 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, 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, + PTRAJPTR& ptraj, KKSTRAWHITCOL& strawhits, KKCALOHITCOL& calohits, KKSTRAWXINGCOL& exings, DOMAINCOL& domains) const; SensorLine 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 @@ -95,8 +107,8 @@ namespace mu2e { StrawResponse const& strawresponse, KKStrawMaterial const& smat, ComboHitCollection const& chcol, Calorimeter const& calo, CCHandle const& cchandle, KKTRK& kktrk) const; - // extend the fit to the surfaces specified in the config - void extendFit(KKTRK& kktrk); + // sample the fit at the specificed surfaces + void sampleFit(KKTRK& kktrk) const; // save the complete fit trajectory as a seed KalSeed createSeed(KKTRK const& kktrk, TrkFitFlag const& seedflag, Calorimeter const& calo, Tracker const& nominalTracker) const; TimeRange range(KKSTRAWHITCOL const& strawhits, KKCALOHITCOL const& calohits, KKSTRAWXINGCOL const& strawxings) const; // time range from a set of hits and element Xings @@ -108,12 +120,11 @@ namespace mu2e { private: void fillTrackerInfo(Tracker const& tracker) const; void addStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf, KKStrawMaterial const& smat, - KKTRK const& kktrk, ComboHitCollection const& chcol, KKSTRAWHITCOL& hits) const; - void addStraws(Tracker const& tracker, KKStrawMaterial const& smat, KKTRK const& kktrk, KKSTRAWHITCOL const& addhits, KKSTRAWXINGCOL& addexings) const; + 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 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 - void extendFit(KKTRK& kktrk) const; - int printLevel_; + int print_; unsigned minNStrawHits_; bool matcorr_, addhits_, addmat_, usecalo_; // flags KKSTRAWHITCLUSTERER shclusterer_; // functor to cluster KKStrawHits @@ -131,22 +142,26 @@ namespace mu2e { // parameters controlling adding hits float maxStrawHitDoca_, maxStrawHitDt_, maxStrawDoca_, maxStrawDocaCon_, maxStrawUposBuff_; int maxDStraw_; // maximum distance from the track a strawhit can be to consider it for adding. - bool skipStrawCheck_; - double addStrawMinDz_; - int strawNBuffer_; // cached info computed from the tracker, used in hit adding; these must be lazy-evaluated as the tracker doesn't exist on construction mutable double strawradius_; mutable double ymin_, ymax_, umax_; // panel-level info mutable double rmin_, rmax_; // plane-level info mutable double spitch_; mutable bool needstrackerinfo_ = true; - - bool saveHitCalib_; + // extrapolation and sampling options + SurfaceMap::SurfacePairCollection sample_; // surfaces to sample the fit + double intertol_; // surface intersection tolerance (mm) + bool sampleinrange_, sampleinbounds_; // require samples to be in range or on surface SaveTraj savetraj_; // trajectory saving option + bool savedomains_; // save domain bounds + bool skipStrawCheck_; + double addStrawMinDz_; + int strawNBuffer_; + bool saveHitCalib_; }; template KKFit::KKFit(KKFitConfig const& fitconfig) : - printLevel_(fitconfig.printLevel()), + print_(fitconfig.printLevel()), minNStrawHits_(fitconfig.minNStrawHits()), matcorr_(fitconfig.matCorr()), addhits_(fitconfig.addHits()), @@ -169,6 +184,10 @@ namespace mu2e { maxStrawDocaCon_(fitconfig.maxStrawDOCAConsistency()), maxStrawUposBuff_(fitconfig.maxStrawUposBuff()), maxDStraw_(fitconfig.maxDStraw()), + intertol_(fitconfig.interTol()), + sampleinrange_(fitconfig.sampleInRange()), + sampleinbounds_(fitconfig.sampleInBounds()), + savedomains_(fitconfig.saveDomains()), skipStrawCheck_(fitconfig.skipStrawCheck()), addStrawMinDz_(fitconfig.addStrawMinDz()), strawNBuffer_(fitconfig.strawNBuffer()), @@ -185,6 +204,15 @@ namespace mu2e { } else { throw cet::exception("RECO")<<"mu2e::KKFit: unknown trajectory option "<< fitconfig.saveTraj() << endl; } + // Lookup surfaces to sample: these should be replaced by extrapolation TODO + SurfaceIdCollection ssids; + for(auto const& sidname : fitconfig.sampleSurfaces()){ + ssids.push_back(SurfaceId(sidname,-1)); // match all elements + } + // translate the sample and extend surface names to actual surfaces using the SurfaceMap. This should come from the + // geometry service eventually, TODO + SurfaceMap smap; + smap.surfaces(ssids,sample_); } template bool KKFit::makeStrawHits(Tracker const& tracker,StrawResponse const& strawresponse,BFieldMap const& kkbf, KKStrawMaterial const& smat, @@ -213,12 +241,81 @@ namespace mu2e { // create the hit. Note these may initially be unusable hits.push_back(std::make_shared(kkbf, pca, combohit, straw, strawidx, strawresponse)); // create the material crossing, including this reference - if(matcorr_) exings.push_back(std::make_shared(hits.back(),smat)); + if(matcorr_){ + auto sline = Mu2eKinKal::strawLine(straw,pca.particleToca()); // line down the straw axis center + CAHint hint(pca.particleToca(),pca.sensorToca()); + PCA spca(ptraj, sline, hint, tprec_ ); + exings.push_back(std::make_shared(hits.back(),spca.localClosestApproach(),smat,straw)); + } if(hits.back()->hitState().usable())ngood++; } return ngood >= minNStrawHits_; } + 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 + PTRAJPTR& ptraj, KKSTRAWHITCOL& strawhits, KKCALOHITCOL& calohits, KKSTRAWXINGCOL& exings, DOMAINCOL& domains) const { // return values + unsigned ngood(0), nactive(0), nsactive(0); + // loop over the TrkStrawHitSeeds in this KalSeed + for(auto const& tshs : kseed.hits()) { + const Straw& straw = tracker.getStraw(tshs.strawId()); + // find the corresponding ComboHit using the index map ( the tshs index is WRT the original ComboHit collection) + auto ifind = strawindexmap.map().find(tshs.index()); + if(ifind == strawindexmap.map().end())throw cet::exception("RECO")<<"mu2e::KKFit: map index not found"<< tshs.index() << endl; + auto chindex = ifind->second; + auto const& combohit = chcol.at(chindex); + auto wline = Mu2eKinKal::hitLine(combohit,straw,strawresponse); // points from the signal to the straw center + // use the recorded TOCA to initialize the POCA + CAHint hint(tshs.particleToca(),tshs.sensorToca()); + PCA pca(*ptraj, wline, hint, tprec_ ); + // create the hit. Note these may initially be unusable + strawhits.push_back(std::make_shared(kkbf, pca, combohit, straw, chindex, strawresponse)); // use original index, so the KSeed can be re-persisted + // set the hit state according to what it was in the fit + strawhits.back()->setState(tshs.wireHitState()); + if(strawhits.back()->hitState().usable())ngood++; + if(strawhits.back()->hitState().active())nactive++; + // create the straw Xing for the associated straw, including the hit reference + } + if(kseed.caloCluster()){ + makeCaloHit(kseed.caloCluster(),calo,*ptraj,calohits); + } + if(matcorr_){ + // add Straw Xings for straws without hits + for(auto const& sx : kseed.straws()){ + auto const& straw = tracker.getStraw(sx._straw); + auto sline = Mu2eKinKal::strawLine(straw,sx._toca); // line down the straw axis center + CAHint hint(sx._toca,sx._toca); + PCA pca(*ptraj, sline, hint, tprec_ ); + // find the associated hit (if any) + KKSTRAWHITPTR shptr; + if(sx.hasHit()){ + for (auto const& sh : strawhits){ + if(sh->strawId() == sx._straw && fabs(sh->time()-pca.particleToca()) < 10.0){// time cude is crude FIXME! + shptr = sh; + break; + } + } + } + exings.push_back(std::make_shared(shptr,pca.localClosestApproach(),smat,straw,sx.active())); + if(sx.active())nsactive++; + } + } + // create domains from the domain boundaries. Use the traj bnom, as that's guaranteed to be consistent + if(kseed.domainBounds().size() > 0){ + for(size_t idb = 0; idb < kseed.domainBounds().size()-1;++idb){ + double tstart = kseed.domainBounds()[idb]; + double trange = kseed.domainBounds()[idb+1] - tstart; + double tmid = tstart + 0.5*trange; + auto domainptr = std::make_shared(tstart,trange,ptraj->nearestPiece(tmid).bnom()); + domains.emplace(domainptr); + } + } + + return ngood >= minNStrawHits_; + } + template SensorLine KKFit::caloAxis(CaloCluster const& cluster, Calorimeter const& calo) const { // move cluster COG into the tracker frame. COG is at the front face of the disk CLHEP::Hep3Vector cog = calo.geomUtil().mu2eToTracker(calo.geomUtil().diskFFToMu2e( cluster.diskID(), cluster.cog3Vector())); @@ -231,7 +328,6 @@ namespace mu2e { return SensorLine(sipmcog,ffcog,cluster.time()+caloDt_,caloPropSpeed_); } - template bool KKFit::makeCaloHit(CCPtr const& cluster, Calorimeter const& calo, PTRAJ const& pktraj, KKCALOHITCOL& hits) const { bool retval(false); auto caxis = caloAxis(*cluster,calo); @@ -262,10 +358,10 @@ namespace mu2e { KKSTRAWHITCOL addstrawhits; KKCALOHITCOL addcalohits; KKSTRAWXINGCOL addstrawxings; - if(addhits_)addStrawHits(tracker, strawresponse, kkbf, smat, kktrk, chcol, addstrawhits ); - if(matcorr_ && addmat_)addStraws(tracker, smat, kktrk, addstrawhits, 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(printLevel_ > 1){ + if(print_ > 1){ std::cout << "KKTrk extension adding " << addstrawhits.size() << " StrawHits and " << addcalohits.size() << " CaloHits and " @@ -275,16 +371,18 @@ namespace mu2e { } template void KKFit::addStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf, KKStrawMaterial const& smat, - KKTRK const& kktrk, ComboHitCollection const& chcol, KKSTRAWHITCOL& addhits) const { - auto const& ftraj = kktrk.fitTraj(); + KKTRK const& kktrk, ComboHitCollection const& chcol, KKSTRAWHITCOL& addhits, KKSTRAWXINGCOL& addexings) const { + auto const& ptraj = kktrk.fitTraj(); // build the set of existing hits std::set oldhits; + std::set oldstrawxs; for(auto const& strawhit : kktrk.strawHits())oldhits.insert(strawhit->strawHitIndex()); + for(auto const& strawx : kktrk.strawXings())oldstrawxs.insert(strawx->strawId()); for( size_t ich=0; ich < chcol.size();++ich){ if(oldhits.find(ich)==oldhits.end()){ // make sure this hit wasn't already found ComboHit const& strawhit = chcol[ich]; if(strawhit.flag().hasAllProperties(addsel_) && (!strawhit.flag().hasAnyProperty(addrej_))){ - double zt = Mu2eKinKal::zTime(ftraj,strawhit.pos().Z(),strawhit.correctedTime()); + double zt = Mu2eKinKal::zTime(ptraj,strawhit.pos().Z(),strawhit.correctedTime()); if(fabs(strawhit.correctedTime()-zt) < maxStrawHitDt_) { // compare the measured time with the estimate from the fit const Straw& straw = tracker.getStraw(strawhit.strawId()); auto wline = Mu2eKinKal::hitLine(strawhit,straw,strawresponse); @@ -292,9 +390,19 @@ namespace mu2e { double htime = wline.measurementTime() - (straw.halfLength()-psign*strawhit.wireDist())/wline.speed(wline.timeAtMidpoint()); CAHint hint(zt,htime); // compute PCA between the trajectory and this straw - PCA pca(ftraj, wline, hint, tprec_ ); + PCA pca(ptraj, wline, hint, tprec_ ); if(fabs(pca.doca()) < maxStrawHitDoca_){ // add test of chi TODO addhits.push_back(std::make_shared(kkbf, pca, strawhit, straw, ich, strawresponse)); + oldhits.insert(addhits.back()->strawHitIndex()); + // add the associated straw + if(matcorr_){ + // test + if(oldstrawxs.find(strawhit.strawId()) != oldstrawxs.end())throw cet::exception("RECO")<<"mu2e::addStrawHits: duplicate straw!!" << std::endl; + auto sline = Mu2eKinKal::strawLine(straw,pca.particleToca()); // line down the straw axis center + CAHint hint(pca.particleToca(),pca.sensorToca()); + PCA spca(ptraj, sline, hint, tprec_ ); + addexings.push_back(std::make_shared(addhits.back(),spca.localClosestApproach(),smat,straw)); + } } } } @@ -303,27 +411,23 @@ namespace mu2e { } template void KKFit::addStraws(Tracker const& tracker, KKStrawMaterial const& smat, KKTRK const& kktrk, - KKSTRAWHITCOL const& addhits, KKSTRAWXINGCOL& addexings) const { + KKSTRAWXINGCOL& addexings) const { // 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 - auto const& ftraj = kktrk.fitTraj(); + auto const& ptraj = kktrk.fitTraj(); // pre-compute some tracker info if needed if(needstrackerinfo_)fillTrackerInfo(tracker); // list the IDs of existing straws: this speeds the search std::set oldstraws; for(auto const& strawxing : kktrk.strawXings())oldstraws.insert(strawxing->strawId()); - // create materials for straw hits just added - for(auto const& strawhit : addhits){ - if(oldstraws.find(strawhit->strawId()) == oldstraws.end()){ - addexings.push_back(std::make_shared(strawhit,smat)); - oldstraws.insert(strawhit->strawId()); - } - } + for(auto const& strawxing : addexings)oldstraws.insert(strawxing->strawId()); // check for vertical tracks - auto t0pos = ftraj.position3(ftraj.t0()); - auto t0dir = ftraj.direction(ftraj.t0()); + auto t0pos = ptraj.position3(ptraj.t0()); + auto t0dir = ptraj.direction(ptraj.t0()); + KKSTRAWHITPTR shptr; // empty ptr if (fabs(t0dir.Z()) < addStrawMinDz_){ + // looping over every plane for vertical tracks is inefficient FIXME for(auto const& plane : tracker.planes()){ if(tracker.planeExists(plane.id())) { for(auto panel_p : plane.panels()){ @@ -332,8 +436,8 @@ namespace mu2e { if (panel.origin().z() > t0pos.Z()+zbuffer || panel.origin().z() < t0pos.Z()-zbuffer){ continue; } - double t = ftraj.t0(); - // for now add all straws + double t = ptraj.t0(); + // test all straws for(unsigned istr = 0; istr < panel.nStraws(); ++istr){ auto const& straw = panel.getStraw(istr); // add strawExists test TODO @@ -342,14 +446,14 @@ namespace mu2e { auto sline = Mu2eKinKal::strawLine(straw,t); // line down the straw axis center CAHint hint(t,t); // compute PCA between the trajectory and this straw - PCA pca(ftraj, sline, hint, tprec_ ); + PCA pca(ptraj, sline, hint, tprec_ ); t = pca.particleToca(); // require consistency with this track passing through this straw double du = fabs((pca.sensorPoca().Vect()-VEC3(straw.wirePosition(0.0))).Dot(VEC3(straw.wireDirection(0.0)))); double doca = fabs(pca.doca()); double dsig = std::max(0.0,doca-strawradius_)/sqrt(pca.docaVar()); if(doca < maxStrawDoca_ && dsig < maxStrawDocaCon_ && du < straw.halfLength() + maxStrawUposBuff_){ - addexings.push_back(std::make_shared(pca.localClosestApproach(),smat,straw)); + addexings.push_back(std::make_shared(shptr,pca.localClosestApproach(),smat,straw)); oldstraws.insert(straw.id()); } } // not existing straw cut @@ -363,8 +467,8 @@ namespace mu2e { if(tracker.planeExists(plane.id())) { double plz = plane.origin().z(); // find the track position in at this plane's central z. - double zt = Mu2eKinKal::zTime(ftraj,plz,ftraj.range().begin()); - auto plpos = ftraj.position3(zt); + double zt = Mu2eKinKal::zTime(ptraj,plz,ptraj.range().begin()); + auto plpos = ptraj.position3(zt); // rough check on the point radius double rho = plpos.Rho(); if((rho > rmin_ && rho < rmax_) || skipStrawCheck_){ @@ -372,7 +476,7 @@ namespace mu2e { for(auto panel_p : plane.panels()){ auto const& panel = *panel_p; // linearly correct position for the track direction due to difference in panel-plane Z position - auto tdir = ftraj.direction(zt); + auto tdir = ptraj.direction(zt); double dz = panel.origin().z()-plz; auto papos = plpos + (dz/tdir.Z())*tdir; // convert this position into panel coordinates @@ -400,13 +504,13 @@ namespace mu2e { auto sline = Mu2eKinKal::strawLine(straw,zt); // line down the straw axis center CAHint hint(zt,zt); // compute PCA between the trajectory and this straw - PCA pca(ftraj, sline, hint, tprec_ ); + PCA pca(ptraj, sline, hint, tprec_ ); // require consistency with this track passing through this straw double du = fabs((pca.sensorPoca().Vect()-VEC3(straw.wirePosition(0.0))).Dot(VEC3(straw.wireDirection(0.0)))); double doca = fabs(pca.doca()); double dsig = std::max(0.0,doca-strawradius_)/sqrt(pca.docaVar()); if(doca < maxStrawDoca_ && dsig < maxStrawDocaCon_ && du < straw.halfLength() + maxStrawUposBuff_){ - addexings.push_back(std::make_shared(pca.localClosestApproach(),smat,straw)); + addexings.push_back(std::make_shared(shptr,pca.localClosestApproach(),smat,straw)); oldstraws.insert(straw.id()); } } // not existing straw cut @@ -420,9 +524,10 @@ namespace mu2e { } // vertical track check } // end function + template void KKFit::addCaloHit(Calorimeter const& calo, KKTRK& kktrk, CCHandle cchandle, KKCALOHITCOL& hits) const { double crystalLength = calo.caloInfo().getDouble("crystalZLength"); - auto const& ftraj = kktrk.fitTraj(); + auto const& ptraj = kktrk.fitTraj(); auto cccol = cchandle.product(); double edep(-1.0); std::shared_ptr chitptr; @@ -435,8 +540,8 @@ namespace mu2e { // test at both faces; if the track is in the right area, test the clusters on this disk // Replace this with an intersection with the calo face TODO for(int iface=0; iface<2; ++iface){ - double zt = Mu2eKinKal::zTime(ftraj,ffpos.z()+iface*crystalLength,ftraj.range().end()); - auto tpos = ftraj.position3(zt); + double zt = Mu2eKinKal::zTime(ptraj,ffpos.z()+iface*crystalLength,ptraj.range().end()); + auto tpos = ptraj.position3(zt); double rho = tpos.Rho(); test[idisk] |= rho > rmin && rho < rmax; } @@ -450,10 +555,10 @@ namespace mu2e { auto caxis = caloAxis(cc,calo); // find the time the seed traj passes the middle of the crystal to form the hint auto pmid = caxis.position3(caxis.timeAtMidpoint()); - double zt = Mu2eKinKal::zTime(ftraj,pmid.Z(),ftraj.range().end()); + double zt = Mu2eKinKal::zTime(ptraj,pmid.Z(),ptraj.range().end()); CAHint hint( zt, caxis.timeAtMidpoint()); // compute closest approach between the fit trajectory and the cluster axis - auto pca = PCA(ftraj, caxis, hint, tprec_ ); + auto pca = PCA(ptraj, caxis, hint, tprec_ ); if(pca.usable() && fabs(pca.doca()) < maxCaloDoca_ && fabs(pca.deltaT()) < maxCaloDt_){ // check that the position is within the active position of the crystal double dz = pca.sensorPoca().Z() - caxis.position3(caxis.measurementTime()).Z(); @@ -512,17 +617,16 @@ namespace mu2e { return TimeRange(tmin,tmax); } - template void KKFit::extendFit(KKTRK& kktrk) { - // extend the fit upstream and downstream (up and down for cosmics) to the specified surfaces; - //TODO - - } - template KalSeed KKFit::createSeed(KKTRK const& kktrk, TrkFitFlag const& seedflag, Calorimeter const& calo, Tracker const& nominalTracker) const { TrkFitFlag fflag(seedflag); // initialize the flag with the seed fit flag if(kktrk.fitStatus().usable()){ fflag.merge(TrkFitFlag::kalmanOK); fflag.merge(TrkFitFlag::seedOK); + fflag.merge(TrkFitFlag::FitOK); + } else { + fflag.clear(TrkFitFlag::kalmanOK); + fflag.clear(TrkFitFlag::seedOK); + fflag.clear(TrkFitFlag::FitOK); } if(kktrk.fitStatus().status_ == Status::converged) fflag.merge(TrkFitFlag::kalmanConverged); if(matcorr_)fflag.merge(TrkFitFlag::MatCorr); @@ -549,19 +653,9 @@ namespace mu2e { if (saveHitCalib_) kseed._hitcalibs.reserve(kktrk.strawHits().size()); for(auto const& strawhit : kktrk.strawHits() ) { - Residual utres = strawhit->refResidual(Mu2eKinKal::tresid); - Residual udres = strawhit->refResidual(Mu2eKinKal::dresid); - Residual ulres = strawhit->refResidual(Mu2eKinKal::lresid); - // compute unbiased residuals; this can fail if the track has marginal coverage - if(kktrk.fitStatus().usable()) { - try { - udres = strawhit->residual(Mu2eKinKal::dresid); - utres = strawhit->residual(Mu2eKinKal::tresid); - ulres = strawhit->residual(Mu2eKinKal::lresid); - } catch (std::exception const& error) { - // std::cout << "Unbiased KKStrawHit residual calculation failure, nDOF = " << fstatus.chisq_.nDOF() << std::endl; - } - } + auto udres = strawhit->residual(Mu2eKinKal::dresid); + auto utres = strawhit->residual(Mu2eKinKal::tresid); + auto ulres = strawhit->residual(Mu2eKinKal::lresid); kseed._hits.emplace_back(strawhit->strawHitIndex(),strawhit->hit(), strawhit->closestApproach().tpData(), strawhit->unbiasedClosestApproach().tpData(), @@ -589,14 +683,7 @@ namespace mu2e { hflag.merge(StrawHitFlag::doca); } // calculate the unbiased time residual - Residual ctres = calohit->refResidual(0); - if(kktrk.fitStatus().usable()){ - try { - ctres = calohit->residual(0); - } catch (std::exception const& error) { - // std::cout << "Unbiased KKCaloHit residual calculation failure, nDOF = " << fstatus.chisq_.nDOF() << std::endl; - } - } + Residual ctres = calohit->residual(0); // calculate the cluster depth = distance along the crystal axis from the POCA to the back face of this disk (where the SiPM sits) double backz = calo.geomUtil().mu2eToTracker(calo.disk(calohit->caloCluster()->diskID()).geomInfo().backFaceCenter()).z(); // calculate the distance from POCA to the SiPM, along the crystal (Z) direction, and projected along the track @@ -633,11 +720,20 @@ namespace mu2e { } // save the fit segments as requested if(kktrk.fitStatus().usable()){ + static const double minrange(1.0e-2); // minimum time range to save a segment. if (savetraj_ == full){ kseed._segments.reserve(fittraj.pieces().size()); for (auto const& traj : fittraj.pieces() ){ - // skip zero-range segments. By convention, sample the state at the mid-time - if(traj->range().range() > 0.0) kseed._segments.emplace_back(*traj,traj->range().mid()); + // skip 'zero-range' segments. By convention, sample the state at the mid-time + if(traj->range().range() > minrange) kseed._segments.emplace_back(*traj,traj->range().mid()); + } + if(savedomains_){ + kseed._domainbounds.reserve(kktrk.domains().size()+1); + for (auto const& domain : kktrk.domains()){ + kseed._domainbounds.push_back(domain->begin()); + } + // save end of last domain + kseed._domainbounds.push_back((*(kktrk.domains().rbegin()))->end()); } } else if (savetraj_ == detector ) { // only save segments inside the tracker volume. First, find the time limits for that @@ -654,7 +750,7 @@ namespace mu2e { tmax = std::max(tmax,inter.time_); } } - // extend as needed to the calohit. Eventually we should also extend to the calorimeter, but that needs to be an extrapolation + // extend as needed to the calohit. Eventually we should extrapolate all tracks to the calorimeter if(kktrk.caloHits().size() > 0){ auto const& calohit = kktrk.caloHits().front(); if(calohit->active()){ @@ -670,10 +766,24 @@ namespace mu2e { kseed._segments.reserve(fittraj.pieces().size());// this will be oversized for (auto const& traj : fittraj.pieces() ){ // skip segments outside the tracker volume range - if(traj->range().range() > 0.0 && (traj->range().inRange(tmin) || traj->range().inRange(tmax) || (traj->range().begin() > tmin && traj->range().end() < tmax)) ) kseed._segments.emplace_back(*traj,traj->range().mid()); + if(traj->range().range() > 0.0 && ((traj->range().begin() > tmin && traj->range().end() < tmax) || traj->range().inRange(tmin) || traj->range().inRange(tmax) ) ) kseed._segments.emplace_back(*traj,traj->range().mid()); } + kseed._segments.shrink_to_fit(); + if(savedomains_){ + kseed._domainbounds.reserve(kktrk.domains().size()+1); // oversized + for (auto const& domain : kktrk.domains()){ + if((domain->range().begin() > tmin && domain->range().end() < tmax) || domain->range().inRange(tmin) || domain->range().inRange(tmax) ) + kseed._domainbounds.push_back(domain->begin()); + // save end of last domain + if(domain->range().inRange(tmax))kseed._domainbounds.push_back(domain->end()); + } + kseed._domainbounds.shrink_to_fit(); + kseed._dtol = kktrk.config().tol_; + } + } else if (savetraj_ == t0seg ) { kseed._segments.emplace_back(t0piece,t0val); + // domain storage is non-sensical in this case, leave empty } } else { kseed._segments.emplace_back(t0piece,t0val); // save the t0 piece even for failed fits @@ -688,24 +798,24 @@ namespace mu2e { template void KKFit::sampleFit(KKTRK const& kktrk,KalIntersectionCollection& inters) const { // add IPA and ST Xings. Sample just before the transit to include the effect of the material static const double epsilon(1.0e-6); - auto const& ftraj = kktrk.fitTraj(); + auto const& ptraj = kktrk.fitTraj(); for(auto const& ipaxing : kktrk.IPAXings()){ double stime = ipaxing->time() - epsilon; - auto const& ktraj = ftraj.nearestPiece(stime); + auto const& ktraj = ptraj.nearestPiece(stime); double dmom,paramomvar,perpmomvar; ipaxing->materialEffects(dmom,paramomvar,perpmomvar); inters.emplace_back(ktraj.stateEstimate(ipaxing->time()),XYZVectorF(ktraj.bnom()),ipaxing->surfaceId(),ipaxing->intersection(),dmom); } for(auto const& stxing : kktrk.STXings()){ double stime = stxing->time() - epsilon; - auto const& ktraj = ftraj.nearestPiece(stime); + auto const& ktraj = ptraj.nearestPiece(stime); double dmom,paramomvar,perpmomvar; stxing->materialEffects(dmom,paramomvar,perpmomvar); inters.emplace_back(ktraj.stateEstimate(stxing->time()),XYZVectorF(ktraj.bnom()),stxing->surfaceId(),stxing->intersection(),dmom); } for(auto const& crvxing : kktrk.CRVXings()){ double stime = crvxing->time() - epsilon; - auto const& ktraj = ftraj.nearestPiece(stime); + auto const& ktraj = ptraj.nearestPiece(stime); double dmom,paramomvar,perpmomvar; crvxing->materialEffects(dmom,paramomvar,perpmomvar); inters.emplace_back(ktraj.stateEstimate(crvxing->time()),XYZVectorF(ktraj.bnom()),crvxing->surfaceId(),crvxing->intersection(),dmom); @@ -714,10 +824,53 @@ namespace mu2e { for(auto const& interpair : kktrk.intersections()) { auto const& sid = std::get<0>(interpair); auto const& inter = std::get<1>(interpair); - auto const& ktraj = ftraj.nearestPiece(inter.time_); + auto const& ktraj = ptraj.nearestPiece(inter.time_); inters.emplace_back(ktraj.stateEstimate(inter.time_),XYZVectorF(ktraj.bnom()),sid,inter); } // sort by time TODO } + + template void KKFit::sampleFit(KKTRK& kktrk) const { + auto const& ptraj = kktrk.fitTraj(); + std::vector ranges; + // test for reflection, and if present, split the test in 2 + auto refltraj = ptraj.reflection(ptraj.range().begin()); + if(refltraj){ + double tmid = refltraj->range().begin(); + ranges.push_back(TimeRange(ptraj.range().begin(),tmid)); + ranges.push_back(TimeRange(tmid,ptraj.range().end())); + } else { + ranges.push_back(ptraj.range()); + } + for(auto range : ranges) { + double tbeg = range.begin(); + double tend = range.end(); + for(auto const& surf : sample_){ + // search for intersections with each surface within the specified time range, going forwards in time + bool goodinter(true); + size_t max_inter = 100; // limit the number of intersections + size_t cur_inter = 0; + // loop to find multiple intersections + while(goodinter && tbeg < tend && cur_inter < max_inter){ + cur_inter += 1; + TimeRange irange(tbeg,tend); + auto surfinter = KinKal::intersect(ptraj,*surf.second,irange,intertol_); + goodinter = surfinter.onsurface_ && ( (! sampleinbounds_) || surfinter.inbounds_ ) && ( (!sampleinrange_) || surfinter.inrange_); + if(goodinter) { + // save the intersection information + kktrk.addIntersection(surf.first,surfinter); + // update for the next intersection + // move past existing intersection to avoid repeating + double epsilon = intertol_/ptraj.speed(surfinter.time_); + tbeg = surfinter.time_ + epsilon; + } + if(print_>1) printf(" [KKFit::%s] Found intersection with surface %15s\n", + __func__, surf.first.name().c_str()); + } + } + } + } + + } #endif diff --git a/Mu2eKinKal/inc/KKFitSettings.hh b/Mu2eKinKal/inc/KKFitSettings.hh index 4b1e7e53a1..9a03aee67c 100644 --- a/Mu2eKinKal/inc/KKFitSettings.hh +++ b/Mu2eKinKal/inc/KKFitSettings.hh @@ -19,10 +19,8 @@ #include "Offline/Mu2eKinKal/inc/StrawXingUpdater.hh" namespace mu2e { namespace Mu2eKinKal{ - - using Name = fhicl::Name; - using Comment = fhicl::Comment; - + using fhicl::Name; + using fhicl::Comment; // struct for defining the KinKal Config object and updaters struct KinKalConfig { fhicl::Atom printLevel { Name("PrintLevel"), Comment("Diagnostic printout Level") }; @@ -92,6 +90,11 @@ namespace mu2e { fhicl::Atom saveHitCalib { Name("SaveHitCalibInfo"), Comment("Save hit calib info in KalSeed") }; // extension and sampling fhicl::Atom saveTraj { Name("SaveTrajectory"), Comment("How to save the trajectory in the KalSeed: None, Full, Detector, or T0 (1 segment containing t0)") }; + fhicl::Atom saveDomains { Name("SaveDomains"), Comment("Save the Bfield domain walls in the KalSeed: this will follow the range specified by SaveTrajectory") }; + fhicl::Sequence sampleSurfaces { Name("SampleSurfaces"), Comment("Sample the final fit at these surfaces") }; + fhicl::Atom interTol { Name("IntersectionTolerance"), Comment("Tolerance for surface intersections (mm)") }; + fhicl::Atom sampleInRange { Name("SampleInRange"), Comment("Require sample times to be inside the fit trajectory time range") }; + fhicl::Atom sampleInBounds { Name("SampleInBounds"), Comment("Require sample intersection point be inside surface bounds (within tolerance)") }; }; // struct for configuring a KinKal fit module struct KKModuleConfig { @@ -103,6 +106,17 @@ namespace mu2e { fhicl::Sequence seederrors { Name("SeedErrors"), Comment("Initial value of seed parameter errors (rms, various units)") }; fhicl::Atom saveAll { Name("SaveAllFits"), Comment("Save all fits, whether they suceed or not"),false }; }; + // Extrapolation configuration + struct KKExtrapConfig { + fhicl::Atom Debug { Name("Debug"), Comment("Debug level"), 0 }; + 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 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") }; + }; } } #endif diff --git a/Mu2eKinKal/inc/KKSHFlag.hh b/Mu2eKinKal/inc/KKSHFlag.hh index 45f9ad8fd0..349fd4477b 100644 --- a/Mu2eKinKal/inc/KKSHFlag.hh +++ b/Mu2eKinKal/inc/KKSHFlag.hh @@ -18,7 +18,8 @@ namespace mu2e { nhdrift=8, // if set, use drift radius to set null hit variance (otherwise use straw radius) longval=9, // if set, use time division to constrain longitudinal position annprob=15, // if set, interpret sign ANN probabilistically - added=20 // record if the hit was added (otherwise original) + added=20, // record if the hit was added (otherwise original) + goodudresid=21, goodutresid, goodulresid // status of unbiased residual calculation }; // functions needed for the BitMap template static std::string const& typeName(); diff --git a/Mu2eKinKal/inc/KKShellXing.hh b/Mu2eKinKal/inc/KKShellXing.hh index 22e58be403..4da6162b6f 100644 --- a/Mu2eKinKal/inc/KKShellXing.hh +++ b/Mu2eKinKal/inc/KKShellXing.hh @@ -22,6 +22,15 @@ namespace mu2e { // construct from a surface, material, intersection, and transverse thickness KKShellXing(SURFPTR surface, SurfaceId const& sid, MatEnv::DetMaterial const& mat, KinKal::Intersection inter, KTRAJPTR reftraj, double thickness, double tol); virtual ~KKShellXing() {} + // clone op for reinstantiation + KKShellXing(KKShellXing const& rhs) = default; + std::shared_ptr< KinKal::ElementXing > clone(CloneContext& context) const override{ + auto rv = std::make_shared< KKShellXing >(*this); + //auto ptr = context.get(reftrajptr_); + auto ptr = std::make_shared(*reftrajptr_); + rv->setReferenceTrajectoryPtr(ptr); + return rv; + }; // ElementXing interface void updateReference(PTRAJ const& ptraj) override; void updateState(MetaIterConfig const& config,bool first) override; @@ -31,10 +40,13 @@ namespace mu2e { KTRAJ const& referenceTrajectory() const override { return *reftrajptr_; } std::vectorconst& matXings() const override { return mxings_; } void print(std::ostream& ost=std::cout,int detail=0) const override; + bool active() const override { return mxings_.size() > 0; } // specific accessors auto const& intersection() const { return inter_; } auto const& material() const { return mat_; } auto const& surfaceId() const { return sid_; } + // other accessors + void setReferenceTrajectoryPtr(KTRAJPTR ptr){ reftrajptr_ = ptr; } private: SURFPTR surf_; // surface SurfaceId sid_; // surface Id diff --git a/Mu2eKinKal/inc/KKStrawHit.hh b/Mu2eKinKal/inc/KKStrawHit.hh index 379a924e36..89d492f49e 100644 --- a/Mu2eKinKal/inc/KKStrawHit.hh +++ b/Mu2eKinKal/inc/KKStrawHit.hh @@ -50,6 +50,34 @@ namespace mu2e { using CA = KinKal::ClosestApproach; KKStrawHit(BFieldMap const& bfield, PCA const& pca, ComboHit const& chit, Straw const& straw, StrawHitIndex const& shindex, StrawResponse const& sresponse); + // clone op for reinstantiation + KKStrawHit(KKStrawHit const& rhs): + bfield_(rhs.bfield()), + whstate_(rhs.hitState()), + dVar_(driftVariance()), + dDdT_(driftVelocity()), + wire_(rhs.wire()), + ca_( + rhs.closestApproach().particleTraj(), + wire_, + rhs.closestApproach().hint(), + rhs.closestApproach().precision() + ), + resids_(rhs.refResiduals()), + chit_(rhs.hit()), + shindex_(rhs.strawHitIndex()), + straw_(rhs.straw()), + sresponse_(rhs.strawResponse()){ + /**/ + }; + std::shared_ptr< KinKal::Hit > clone(CloneContext& context) const override{ + auto rv = std::make_shared< KKStrawHit >(*this); + auto ca = rv->closestApproach(); + auto trajectory = std::make_shared(ca.particleTraj()); + ca.setTrajectory(trajectory); + rv->setClosestApproach(ca); + return rv; + }; // Hit interface implementations void updateState(MetaIterConfig const& config,bool first) override; unsigned nResid() const override { return 3; } // potentially 2 residuals @@ -62,6 +90,8 @@ namespace mu2e { void updateReference(PTRAJ const& ptraj) override; KTRAJPTR const& refTrajPtr() const override { return ca_.particleTrajPtr(); } void print(std::ostream& ost=std::cout,int detail=0) const override; + // re-override; even though this is implemented in the base class + bool active() const override { return whstate_.active(); } // accessors auto const& closestApproach() const { return ca_; } auto const& hitState() const { return whstate_; } @@ -79,6 +109,8 @@ namespace mu2e { auto updater() const { return whstate_.algo_; } void setState(WireHitState const& whstate); // allow cluster updaters to set the state directly DriftInfo fillDriftInfo(CA const& ca) const; + auto const& driftVariance() { return dVar_; } + auto const& driftVelocity() { return dDdT_; } private: BFieldMap const& bfield_; // drift calculation requires the BField for ExB effects WireHitState whstate_; // current state @@ -97,6 +129,8 @@ namespace mu2e { StrawResponse const& sresponse_; // straw calibration information // utility functions void updateWHS(MetaIterConfig const& miconfig); + // clone support + void setClosestApproach(const CA& ca){ ca_ = ca; } }; // struct to sort hits by time @@ -189,8 +223,12 @@ namespace mu2e { } template DriftInfo KKStrawHit::fillDriftInfo(CA const& ca) const { - double lorentzAngle = Mu2eKinKal::LorentzAngle(ca.tpData(),ca.particleTraj().bnom().Unit()); - return sresponse_.driftInfo(strawId(),ca.deltaT(),lorentzAngle); + if(ca.usable()){ + double lorentzAngle = Mu2eKinKal::LorentzAngle(ca.tpData(),ca.particleTraj().bnom().Unit()); + return sresponse_.driftInfo(strawId(),ca.deltaT(),lorentzAngle); + } else { + return DriftInfo(); + } } template VEC3 KKStrawHit::dRdX(unsigned ires) const { @@ -215,22 +253,22 @@ namespace mu2e { if(whstate.constrainTOT()){ double tvar = chit_.timeVar(); double dt = ca_.deltaT() - chit_.driftTime(); - resids[Mu2eKinKal::tresid] = Residual(dt,tvar,0.0,true,ca_.dTdP()); + resids[Mu2eKinKal::tresid] = Residual(dt,tvar,0.0,ca_.dTdP()); } // distance residual if(whstate.driftConstraint()){ double dr = whstate.lrSign()*dinfo.rDrift_ - ca_.doca(); DVEC dRdP = whstate.lrSign()*dDdT_*ca_.dTdP() -ca_.dDdP(); - resids[Mu2eKinKal::dresid] = Residual(dr,dVar_,0.0,true,dRdP); + resids[Mu2eKinKal::dresid] = Residual(dr,dVar_,0.0,dRdP); } else { // Null LR ambiguity. interpret DOCA against the wire directly as the spatial residual - resids[Mu2eKinKal::dresid] = Residual(ca_.doca(),dVar_,0.0,true,ca_.dDdP()); + resids[Mu2eKinKal::dresid] = Residual(ca_.doca(),dVar_,0.0,ca_.dDdP()); // optionally use the null hit time measurement to constrain t0 if(whstate.constrainAbsDriftDt()){ double dt = ca_.deltaT() - sresponse_.strawDrift().D2T(fabs(ca_.doca()),dinfo.LorentzAngle_); double tvar = dinfo.driftTimeVar(); // this overwrites the TOT time constraint, in principle both can be used TODO - resids[Mu2eKinKal::tresid] = Residual(dt,tvar,0.0,true,ca_.dTdP()); + resids[Mu2eKinKal::tresid] = Residual(dt,tvar,0.0,ca_.dTdP()); } } @@ -256,7 +294,7 @@ namespace mu2e { double pdir_dot_sdir = ca_.sensorDirection().Dot(ca_.particleDirection()); DVEC dLdP = dx_dot_sdir + (dx_dot_sdir*pdir_dot_sdir - dx_dot_pdir - d_dot_dm)/(1-pdir_dot_sdir*pdir_dot_sdir)*pdir_dot_sdir; dLdP *= -1; - resids[Mu2eKinKal::lresid] = Residual(lresidval,lresidvar,0.0,true,dLdP); + resids[Mu2eKinKal::lresid] = Residual(lresidval,lresidvar,0.0,dLdP); } } diff --git a/Mu2eKinKal/inc/KKStrawHitCluster.hh b/Mu2eKinKal/inc/KKStrawHitCluster.hh index 979690d80b..a3ccbb0047 100644 --- a/Mu2eKinKal/inc/KKStrawHitCluster.hh +++ b/Mu2eKinKal/inc/KKStrawHitCluster.hh @@ -54,6 +54,22 @@ namespace mu2e { KKStrawHitCluster(KKSTRAWHITPTR const& hitptr); // create from a collection of panel hits KKStrawHitCluster(KKSTRAWHITCOL const& hits,KKSTRAWHITCLUSTERER const& clusterer); + // clone op for reinstantiation + KKStrawHitCluster(KKStrawHitCluster const& rhs){ + /**/ + }; + std::shared_ptr< KinKal::Hit > clone(CloneContext& context) const override{ + auto rv = std::make_shared< KKStrawHitCluster >(*this); + for (const auto& ptr: this->strawHits()){ + auto hit = context.get(ptr); + rv->push_back(hit); + } + for (const auto& ptr: this->strawXings()){ + auto xng = context.get(ptr); + rv->push_back(xng); + } + return rv; + }; //Hit interface bool active() const override { return false; } // panel hits are never active KinKal::Chisq chisq(KinKal::Parameters const& params) const override { return KinKal::Chisq(); } @@ -73,6 +89,11 @@ namespace mu2e { bool canAddHit(KKSTRAWHITPTR hit,KKSTRAWHITCLUSTERER const& clusterer) const; void addHit(KKSTRAWHITPTR hit,KKSTRAWHITCLUSTERER const& clusterer); void addXing(KKSTRAWXINGPTR xing); + + protected: + void push_back(KKSTRAWHITPTR hit){ hits_.push_back(hit); } + void push_back(KKSTRAWXINGPTR xng){ xings_.push_back(xng); }; + private: // references to the individual hits and xings in this hit cluster KKSTRAWHITCOL hits_; diff --git a/Mu2eKinKal/inc/KKStrawXing.hh b/Mu2eKinKal/inc/KKStrawXing.hh index dd51c54364..85dc464e74 100644 --- a/Mu2eKinKal/inc/KKStrawXing.hh +++ b/Mu2eKinKal/inc/KKStrawXing.hh @@ -1,7 +1,7 @@ #ifndef Mu2eKinKal_KKStrawXing_hh #define Mu2eKinKal_KKStrawXing_hh // -// StrawXing using Mu2e-specific StrawMaterial class. Otherwise it's the same as KinKal::StrawXing +// StrawXing using Mu2e-specific StrawMaterial class // #include "KinKal/Detector/ElementXing.hh" #include "Offline/Mu2eKinKal/inc/KKStrawHit.hh" @@ -29,11 +29,42 @@ namespace mu2e { using CA = KinKal::ClosestApproach; using KKSTRAWHIT = KKStrawHit; using KKSTRAWHITPTR = std::shared_ptr; - // construct without an associated StrawHit - KKStrawXing(CA const& ca, KKStrawMaterial const& smat, Straw const& straw); + // construct with closest approach + KKStrawXing(KKSTRAWHITPTR const& strawhit, CA const& ca, KKStrawMaterial const& smat, Straw const& straw,bool active=false); // construct with an associated StrawHit KKStrawXing(KKSTRAWHITPTR const& strawhit, KKStrawMaterial const& smat); virtual ~KKStrawXing() {} + // clone op for reinstantiation + KKStrawXing(KKStrawXing const& rhs): + KKStrawXing( + rhs.strawHitPtr(), + rhs.closestApproach(), + rhs.strawMaterial(), + rhs.straw() + ){ + auto shptr = rhs.strawHitPtr(); + if (shptr){ + this->setStrawHitPtr(rhs.strawHitPtr()); + } + } + std::shared_ptr< KinKal::ElementXing > clone(CloneContext& context) const override{ + auto rv = std::make_shared< KKStrawXing >(*this); + + // point to new instance of partner hit, if not null + KKSTRAWHITPTR shptr; + if (shptr_){ + shptr = context.get(shptr_); + } + rv->setStrawHitPtr(shptr); + + // point to new instance of ClosestApproach + auto ca = rv->closestApproach(); + auto trajectory = std::make_shared(ca.particleTraj()); + ca.setTrajectory(trajectory); + rv->setClosestApproach(ca); + + return rv; + }; // ElementXing interface void updateReference(PTRAJ const& ptraj) override; void updateState(MetaIterConfig const& config,bool first) override; @@ -45,17 +76,17 @@ namespace mu2e { double transitTime() const override; // time to cross this element KTRAJ const& referenceTrajectory() const override { return ca_.particleTraj(); } void print(std::ostream& ost=std::cout,int detail=0) const override; + bool active() const override; // accessors - bool active() const { return active_; } + auto const& config() const { return sxconfig_; } auto const& closestApproach() const { return ca_; } auto const& strawMaterial() const { return smat_; } - auto const& config() const { return sxconfig_; } auto precision() const { return ca_.precision(); } auto const& straw() const { return straw_; } auto const& strawId() const { return straw_.id(); } auto const& strawHitPtr() const { return shptr_; } private: - StrawXingUpdater sxconfig_; // cache of most recent config + StrawXingUpdater sxconfig_; // cache of most recent update KKSTRAWHITPTR shptr_; // reference to associated StrawHit SensorLine axis_; // straw axis, expressed as a timeline KKStrawMaterial const& smat_; @@ -64,17 +95,21 @@ namespace mu2e { double toff_; // small time offset std::vector mxings_; Parameters fparams_; // parameter change for forwards time - bool active_; // active or not + bool active_; // inside active region or not double varscale_; // variance scale + // modifiers to support cloning + void setStrawHitPtr(KKSTRAWHITPTR ptr) { shptr_ = ptr; } + void setClosestApproach(const CA& ca){ ca_ = ca; } }; - template KKStrawXing::KKStrawXing(CA const& ca, KKStrawMaterial const& smat, Straw const& straw) : + template KKStrawXing::KKStrawXing(KKSTRAWHITPTR const& strawhit, CA const& ca, KKStrawMaterial const& smat, Straw const& straw,bool active) : + shptr_(strawhit), axis_(ca.sensorTraj()), smat_(smat), straw_(straw), ca_(ca.particleTraj(),axis_,ca.hint(),ca.precision()), toff_(smat.wireRadius()/ca.particleTraj().speed(ca.particleToca())), // locate the effect to 1 side of the wire to avoid overlap with hits - active_(false), + active_(active), varscale_(1.0) {} @@ -89,6 +124,11 @@ namespace mu2e { varscale_(1.0) {} + template bool KKStrawXing::active() const { + // if the associated hit is active, use it's state. Otherwise use the intrinsic state + return active_ || (shptr_ && shptr_->hitState().active()); + } + template void KKStrawXing::updateReference(PTRAJ const& ptraj) { CAHint tphint(axis_.timeAtMidpoint(),axis_.timeAtMidpoint()); if(ca_.usable()){ @@ -106,28 +146,26 @@ namespace mu2e { template void KKStrawXing::updateState(MetaIterConfig const& miconfig,bool first) { // reset - fparams_ = Parameters(); + mxings_.clear(); if(first) { - active_ = false; - // search for an update to the xing configuration among this meta-iteration payload + // search for an updater among this meta-iteration payload auto sxconfig = miconfig.findUpdater(); - if(sxconfig != 0){ - sxconfig_ = *sxconfig; - } + // cache, as this also sets parameters used in calculating path lengths + if(sxconfig != 0)sxconfig_ = *sxconfig; + // require a validat updater + if(sxconfig_.nsig_ < 0)throw cet::exception("RECO")<<"mu2e::KKStrawXing: invalid updater!" << std::endl; + // update the DOCA range test + if(sxconfig_.maxdoca_ > 0.0)active_ = fabs(ca_.tpData().doca()) < sxconfig_.maxdoca_; + // set the variance scale (temperature) if(sxconfig_.scalevar_) varscale_ = miconfig.varianceScale(); else varscale_ = 1.0; - // update the associated hit state - if(shptr_)active_ = shptr_->hitState().active(); - // decide if this straw xing should be active - active_ |= fabs(ca_.tpData().doca()) < sxconfig_.maxdoca_; - } - if(active_){ - // update the material xings from gas, straw wall, and wire - smat_.findXings(ca_.tpData(),sxconfig_,mxings_); - fparams_ = this->parameterChange(varscale_); } + // update the material xings from gas, straw wall, and wire + smat_.findXings(ca_.tpData(),sxconfig_,mxings_); + // update the effect these have on the parameters + fparams_ = this->parameterChange(varscale_); } template double KKStrawXing::transitTime() const { diff --git a/Mu2eKinKal/inc/KKTrack.hh b/Mu2eKinKal/inc/KKTrack.hh index 5ca8a5cefd..c6dc11848c 100644 --- a/Mu2eKinKal/inc/KKTrack.hh +++ b/Mu2eKinKal/inc/KKTrack.hh @@ -54,9 +54,68 @@ namespace mu2e { using KKINTER = std::tuple; using KKINTERCOL = std::vector; using TRACK = KinKal::Track; + using PKTRAJ = KinKal::ParticleTrajectory; + using PKTRAJPTR = std::unique_ptr; + using DOMAINPTR = std::shared_ptr; + using DOMAINCOL = std::set; // construct from configuration, fit environment, and hits and materials KKTrack(Config const& config, BFieldMap const& bfield, KTRAJ const& seedtraj, PDGCode::type tpart, KKSTRAWHITCLUSTERER const& shclusterer, KKSTRAWHITCOL const& strawhits, KKSTRAWXINGCOL const& strawxings, KKCALOHITCOL const& calohits, std::array constraints = {0}); + // construct from regrown constituants + KKTrack(Config const& config, BFieldMap const& bfield, PDGCode::type tpart, + PKTRAJPTR& fittraj, KKSTRAWHITCOL& strawhits, KKSTRAWXINGCOL& strawxings, KKCALOHITCOL& calohits, DOMAINCOL& domains); + // copy constructor + KKTrack(KKTrack const& rhs, CloneContext& context): KinKal::Track(rhs, context), + tpart_(rhs.fitParticle()), + shclusterer_(rhs.strawHitClusterer()), + strawhits_(rhs.strawHits()), + strawxings_(rhs.strawXings()), + inters_(rhs.intersections()), + calohits_(rhs.caloHits()){ + // hits and crossings were reallocated into the base copy; here, + // we propagate references to those reallocations to the subclass + this->remap_pointer_collection(strawhits_, this->hits(), context); + this->remap_pointer_collection(calohits_, this->hits(), context); + this->remap_pointer_collection(strawxings_, this->exings(), context); + // noncontained crossings and clusters are not stored in + // the base, so we are responsible for reallocations here + this->clone_pointer_collection(ipaxings_, rhs.IPAXings(), context); + this->clone_pointer_collection(stxings_, rhs.STXings(), context); + this->clone_pointer_collection(crvxings_, rhs.CRVXings(), context); + this->clone_pointer_collection(strawhitclusters_, rhs.strawHitClusters(), context); + } + + template + using PtrVector = std::vector< std::shared_ptr >; + template + void remap_pointer_collection(PtrVector& lhs, const PtrVector& rhs, + CloneContext& context){ + lhs.clear(); + lhs.reserve(rhs.size()); + for (const auto& item: rhs){ + auto realloc = context.get(item); + auto casted = dynamic_pointer_cast(realloc); + if (casted){ + lhs.push_back(casted); + } + } + } + template + void clone_pointer_collection(PtrVector& lhs, const PtrVector& rhs, + CloneContext& context){ + lhs.clear(); + lhs.reserve(rhs.size()); + for (const auto& ptr: rhs){ + auto realloc = ptr->clone(context); + auto casted = dynamic_pointer_cast(realloc); + if (!casted){ + std::string msg = "mistyped pointee"; + throw cet::exception("KKTrack") << msg << std::endl; + } + lhs.push_back(casted); + } + } + // extend the track according to new configuration, hits, and/or exings void extendTrack(Config const& config, KKSTRAWHITCOL const& strawhits, KKSTRAWXINGCOL const& strawxings, KKCALOHITCOL const& calohits ); @@ -84,18 +143,21 @@ namespace mu2e { void printFit(std::ostream& ost=std::cout,int detail=0) const; // allow reversing the charge (in some fits it can change dynamically) void reverseCharge() { tpart_ = static_cast(tpart_*-1); } + // other accessors + auto const& strawHitClusterer() const { return shclusterer_; } private: // record the particle type PDGCode::type tpart_; KKSTRAWHITCLUSTERER shclusterer_; // clustering functor + TrkFitFlag flag_; KKSTRAWHITCOL strawhits_; // straw hits used in this fit KKSTRAWXINGCOL strawxings_; // straw material crossings used in this fit + KKCALOHITCOL calohits_; // calo hits used in this fit KKIPAXINGCOL ipaxings_; // ipa material crossings used in extrapolation KKSTXINGCOL stxings_; // stopping target material crossings used in extrapolation KKCRVXINGCOL crvxings_; // crv crossings using in extrapolation KKINTERCOL inters_; // other recorded intersections KKSTRAWHITCLUSTERCOL strawhitclusters_; // straw hit clusters used in this fit - KKCALOHITCOL calohits_; // calo hits used in this fit // utility function to convert to generic types void convertTypes( KKSTRAWHITCOL const& strawhits, KKSTRAWXINGCOL const& strawxings,KKCALOHITCOL const& calohits, MEASCOL& hits, EXINGCOL& exings); @@ -109,10 +171,8 @@ namespace mu2e { KKSTRAWXINGCOL const& strawxings, KKCALOHITCOL const& calohits, std::array constraints) : - KinKal::Track(config,bfield,seedtraj), tpart_(tpart), shclusterer_(shclusterer), - strawhits_(strawhits), - strawxings_(strawxings), - calohits_(calohits) { + KinKal::Track(config,bfield), tpart_(tpart), shclusterer_(shclusterer), + strawhits_(strawhits), strawxings_(strawxings), calohits_(calohits) { MEASCOL hits; // polymorphic container of hits EXINGCOL exings; // polymorphic container of detector element crossings // add the hits to clusters, as required @@ -123,7 +183,7 @@ namespace mu2e { shcluster->print(std::cout,this->config().plevel_); } } - convertTypes(strawhits_, strawxings_, calohits_, hits,exings); + convertTypes(strawhits_, strawxings_, calohits_, hits, exings); std::array mask = {false}; bool constraining = false; @@ -148,10 +208,23 @@ namespace mu2e { } hits.push_back(std::make_shared>(seedtraj.range().mid(),seedtraj,cparams,mask)); } - - this->fit(hits,exings); + // now fit these + this->fit(hits, exings, seedtraj); } + + template KKTrack::KKTrack(Config const& config, BFieldMap const& bfield, PDGCode::type tpart, + PKTRAJPTR& fittraj, KKSTRAWHITCOL& strawhits, KKSTRAWXINGCOL& strawxings, KKCALOHITCOL& calohits, DOMAINCOL& domains) : + KinKal::Track(config,bfield), tpart_(tpart),shclusterer_(StrawIdMask::none,0,0.0), + strawhits_(strawhits), strawxings_(strawxings), calohits_(calohits) { + // convert types + MEASCOL hits; // polymorphic container of hits + EXINGCOL exings; // polymorphic container of detector element crossings + convertTypes(strawhits_, strawxings_, calohits_, hits, exings); + this->fit(hits,exings,domains,fittraj); + // add ParameterHit, intersect, extrapolate, ... TODO + } + template void KKTrack::addHitClusters(KKSTRAWHITCOL const& strawhits,KKSTRAWXINGCOL const& strawxings,MEASCOL& hits) { if(shclusterer_.clusterLevel() != StrawIdMask::none){ for(auto const& strawhitptr : strawhits){ @@ -238,7 +311,6 @@ namespace mu2e { } template void KKTrack::printFit(std::ostream& ost,int printlevel) const { - if(printlevel > 1) std::cout << "Seed Helix " << this->seedTraj() << std::endl; TRACK::print(ost,0); ost << "Fit with " << strawhits_.size() << " StrawHits and " << calohits_.size() << " CaloHits and " << strawxings_.size() << " Straw Xings" << std::endl; if(printlevel > 2){ @@ -276,7 +348,5 @@ namespace mu2e { crvxings_.push_back(crvxingptr); } - - } #endif diff --git a/Mu2eKinKal/inc/StrawXingUpdater.hh b/Mu2eKinKal/inc/StrawXingUpdater.hh index 4fb737003c..5de8a255df 100644 --- a/Mu2eKinKal/inc/StrawXingUpdater.hh +++ b/Mu2eKinKal/inc/StrawXingUpdater.hh @@ -7,9 +7,9 @@ namespace mu2e { struct StrawXingUpdater { using SXUConfig = std::tuple; static std::string const& configDescription(); // description of the variables - double maxdoca_ =0; // maximum DOCA to include this straw's material - double nsig_ =0; // Number of doca_sigma around doca value to use when averageing - bool scalevar_ =false; // scale variance or not + double maxdoca_ =0; // maximum DOCA to activate straw materials without an associated hits. <0 means don't change the state + double nsig_ =-1; // Number of doca_sigma around doca value to use when averageing + bool scalevar_ =false; // scale variance with temperature or not int diag_ =0; // diag print level // default constructor is functional but will always use the impact-parameter averaged material StrawXingUpdater(SXUConfig const& sxusetting); diff --git a/Mu2eKinKal/src/KKSHFlag.cc b/Mu2eKinKal/src/KKSHFlag.cc index 738a0eabc1..cb9271d1c1 100644 --- a/Mu2eKinKal/src/KKSHFlag.cc +++ b/Mu2eKinKal/src/KKSHFlag.cc @@ -23,6 +23,9 @@ namespace mu2e { bitnames[std::string("LongVal")] = bit_to_mask(longval); bitnames[std::string("ANNProb")] = bit_to_mask(annprob); bitnames[std::string("Added")] = bit_to_mask(added); + bitnames[std::string("GoodUDResid")] = bit_to_mask(goodudresid); + bitnames[std::string("GoodUTResid")] = bit_to_mask(goodutresid); + bitnames[std::string("GoodULResid")] = bit_to_mask(goodulresid); } return bitnames; } diff --git a/Mu2eKinKal/src/KKStrawMaterial.cc b/Mu2eKinKal/src/KKStrawMaterial.cc index cb033bb556..a79a54cef3 100644 --- a/Mu2eKinKal/src/KKStrawMaterial.cc +++ b/Mu2eKinKal/src/KKStrawMaterial.cc @@ -40,7 +40,7 @@ namespace mu2e { wallpath = gaspath = wirepath = 0.0; PathCalc retval = KKStrawMaterial::unknown; double docarange = caconfig.nsig_*sqrt(std::max(0.0,cadata.docaVar())); - // if the doca range covers the straw, use averages + // 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); double imaxdoca = std::min(irad_,adoca+docarange); diff --git a/Mu2eKinKal/src/LoopHelixFit_module.cc b/Mu2eKinKal/src/LoopHelixFit_module.cc index 769a3a8607..e9aec763d1 100644 --- a/Mu2eKinKal/src/LoopHelixFit_module.cc +++ b/Mu2eKinKal/src/LoopHelixFit_module.cc @@ -64,10 +64,7 @@ #include "Offline/Mu2eKinKal/inc/KKBField.hh" #include "Offline/Mu2eKinKal/inc/KKConstantBField.hh" #include "Offline/Mu2eKinKal/inc/KKFitUtilities.hh" -#include "Offline/Mu2eKinKal/inc/ExtrapolateToZ.hh" -#include "Offline/Mu2eKinKal/inc/ExtrapolateIPA.hh" -#include "Offline/Mu2eKinKal/inc/ExtrapolateST.hh" -#include "Offline/Mu2eKinKal/inc/KKShellXing.hh" +#include "Offline/Mu2eKinKal/inc/KKExtrap.hh" // C++ #include #include @@ -75,8 +72,6 @@ #include #include // -// Original author D. Brown (LBNL) 11/18/20 -// namespace mu2e { using KTRAJ= KinKal::LoopHelix; using PTRAJ = KinKal::ParticleTrajectory; @@ -88,9 +83,6 @@ namespace mu2e { using KKSTRAWXING = KKStrawXing; using KKSTRAWXINGPTR = std::shared_ptr; using KKSTRAWXINGCOL = std::vector; - using KKIPAXING = KKShellXing; - using KKIPAXINGPTR = std::shared_ptr; - using KKIPAXINGCOL = std::vector; using KKSTXING = KKShellXing; using KKSTXINGPTR = std::shared_ptr; using KKSTXINGCOL = std::vector; @@ -133,19 +125,6 @@ namespace mu2e { struct KKLHModuleConfig : KKModuleConfig { fhicl::Sequence seedCollections {Name("HelixSeedCollections"), Comment("Seed fit collections to be processed ") }; fhicl::OptionalAtom fixedBField { Name("ConstantBField"), Comment("Constant BField value") }; - fhicl::Sequence sampleSurfaces { Name("SampleSurfaces"), Comment("When creating the KalSeed, sample the fit at these surfaces") }; - fhicl::Atom sampleInRange { Name("SampleInRange"), Comment("Require sample times to be inside the fit trajectory time range") }; - fhicl::Atom sampleInBounds { Name("SampleInBounds"), Comment("Require sample intersection point be inside surface bounds (within tolerance)") }; - fhicl::Atom interTol { Name("IntersectionTolerance"), Comment("Tolerance for surface intersections (mm)") }; - }; - // Extrapolation configuration - struct KKExtrapConfig { - fhicl::Atom MaxDt { Name("MaxDt"), Comment("Maximum time 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 Debug { Name("Debug"), Comment("Debug level"), 0 }; }; struct HelixMaskConfig { fhicl::OptionalAtom minHelixP{ Name("MinHelixMom"), Comment("Minimum Momentum of a helix for a track to be fit.")}; @@ -157,7 +136,7 @@ namespace mu2e { fhicl::Table extSettings { Name("ExtensionSettings") }; fhicl::Table matSettings { Name("MaterialSettings") }; fhicl::OptionalTable finalSettings { Name("FinalSettings") }; - fhicl::OptionalTable Extrapolation { Name("Extrapolation") }; + fhicl::OptionalTable extrapSettings { Name("ExtrapolationSettings") }; // LoopHelix module specific config fhicl::OptionalAtom slopeSigThreshold{ Name("SlopeSigThreshold"), Comment("Input helix seed slope significance threshold to assume the direction")}; fhicl::OptionalAtom fitDirection { Name("FitDirection"), Comment("Particle direction to fit, either \"upstream\" or \"downstream\"")}; @@ -176,19 +155,10 @@ namespace mu2e { TrkFitFlag fitflag_; // parameter-specific functions that need to be overridden in subclasses KTRAJ makeSeedTraj(HelixSeed const& hseed,TimeRange const& trange,VEC3 const& bnom, int charge) const; - bool goodFit(KKTRK const& ktrk) const; + bool goodFit(KKTRK const& ktrk,KTRAJ const& seed) const; bool goodHelix(HelixSeed const& hseed) const; std::vector chooseHelixDir(HelixSeed const& hseed) const; std::unique_ptr fitTrack(art::Event& event, HelixSeed const& hseed, const TrkFitDirection fdir, const PDGCode::type fitpart); - // extrapolation functions - void extrapolate(KKTRK& ktrk) const; - void toTrackerEnds(KKTRK& ktrk) const; - bool extrapolateIPA(KKTRK& ktrk,TimeDir trkdir) const; - bool extrapolateST(KKTRK& ktrk,TimeDir trkdir) const; - bool extrapolateTracker(KKTRK& ktrk,TimeDir tdir) const; - bool extrapolateTSDA(KKTRK& ktrk,TimeDir tdir) const; - void toOPA(KKTRK& ktrk, double tstart, TimeDir tdir) const; - void sampleFit(KKTRK& kktrk) const; void print_track_info(const KalSeed& kkseed, const KKTRK& ktrk) const; // data payload @@ -214,23 +184,10 @@ namespace mu2e { Config config_; // initial fit configuration object Config exconfig_; // extension configuration object Config fconfig_; // final final configuration object - double intertol_; // surface intersection tolerance (mm) - bool sampleinrange_, sampleinbounds_; // require samples to be in range or on surface + std::unique_ptr extrap_; // extrapolation helper bool fixedfield_; // special case usage for seed fits, if no BField corrections are needed - SurfaceMap smap_; - AnnPtr tsdaptr_; - DiskPtr trkfrontptr_, trkmidptr_, trkbackptr_; - FruPtr opaptr_; - bool extrapolate_, backToTracker_, toOPA_, toTrackerEnds_, upstream_; - ExtrapolateToZ TSDA_, trackerFront_, trackerBack_; // extrapolation predicate based on Z values - ExtrapolateIPA extrapIPA_; // extrapolation to intersections with the IPA - ExtrapolateST extrapST_; // extrapolation to intersections with the ST - 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 - SurfaceMap::SurfacePairCollection sample_; // surfaces to sample the fit //Helix Mask params float minHelixP_ = -1.; - CylPtr STInner_, STOuter_; int nSeen_ = 0; int nFit_ = 0; int nSkipped_ = 0; @@ -253,10 +210,7 @@ namespace mu2e { kkmat_(settings().matSettings()), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), exconfig_(Mu2eKinKal::makeConfig(settings().extSettings())), - intertol_(settings().modSettings().interTol()), - sampleinrange_(settings().modSettings().sampleInRange()), - sampleinbounds_(settings().modSettings().sampleInBounds()), - fixedfield_(false), extrapolate_(false), backToTracker_(false), toOPA_(false) + fixedfield_(false) { std::string fdir; if(settings().fitDirection(fdir))fdir_ = fdir; @@ -278,6 +232,8 @@ namespace mu2e { fixedfield_ = true; kkbf_ = std::move(std::make_unique(VEC3(0.0,0.0,bz))); } + // setup extrapolation + if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings(),kkmat_); // setup optional fit finalization; this just updates the internals, not the fit result itself if(settings().finalSettings()){ @@ -290,54 +246,11 @@ namespace mu2e { fconfig_.convdchisq_ = settings().finalSettings()->convdchisq(); fconfig_.maxniter_ = settings().finalSettings()->maxniter(); } - // configure extrapolation - if(settings().Extrapolation()){ - extrapolate_ = true; - backToTracker_ = settings().Extrapolation()->BackToTracker(); - toTrackerEnds_ = settings().Extrapolation()->ToTrackerEnds(); - upstream_ = settings().Extrapolation()->Upstream(); - toOPA_ = settings().Extrapolation()->ToOPA(); - auto const& IPA = smap_.DS().innerProtonAbsorberPtr(); - // global configs - double maxdt = settings().Extrapolation()->MaxDt(); - double btol = settings().extSettings().btol(); // use the same BField cor. tolerance as in fit extension - int debug = settings().Extrapolation()->Debug(); - // predicate to extrapolate through IPA - extrapIPA_ = ExtrapolateIPA(maxdt,btol,intertol_,IPA,debug); - // predicate to extrapolate through ST - if(debug > 0)std::cout << "IPA limits z " << extrapIPA_.zmin() << " " << extrapIPA_.zmax() << std::endl; - extrapST_ = ExtrapolateST(maxdt,btol,intertol_,smap_.ST(),debug); - // temporary - if(debug > 0)std::cout << "ST limits z " << extrapST_.zmin() << " " << extrapST_.zmax() << " r " << extrapST_.rmin() << " " << extrapST_.rmax() << std::endl; - // extrapolate to the front of the tracker - trackerFront_ = ExtrapolateToZ(maxdt,btol,smap_.tracker().front().center().Z(),debug); - trackerBack_ = ExtrapolateToZ(maxdt,btol,smap_.tracker().back().center().Z(),debug); - // extrapolate to the back of the detector solenoid - TSDA_ = ExtrapolateToZ(maxdt,btol,smap_.DS().upstreamAbsorber().center().Z(),debug); - tsdaptr_ = smap_.DS().upstreamAbsorberPtr(); - trkfrontptr_ = smap_.tracker().frontPtr(); - trkmidptr_ = smap_.tracker().middlePtr(); - trkbackptr_ = smap_.tracker().backPtr(); - opaptr_ = smap_.DS().outerProtonAbsorberPtr(); - } if (settings().HelixMask()){ if (settings().HelixMask()->minHelixP()) {minHelixP_ = settings().HelixMask()->minHelixP().value();} } - - // additional surfaces to sample: these should be replaced by extrapolation TODO - SurfaceIdCollection ssids; - for(auto const& sidname : settings().modSettings().sampleSurfaces()){ - ssids.push_back(SurfaceId(sidname,-1)); // match all elements - } - // translate the sample and extend surface names to actual surfaces using the SurfaceMap. This should come from the - // geometry service eventually, TODO - SurfaceMap smap; - smap.surfaces(ssids,sample_); - // Find the target bounding surfaces as well - STInner_ = smap.ST().innerPtr(); - STOuter_ = smap.ST().outerPtr(); } void LoopHelixFit::beginRun(art::Run& run) { @@ -417,7 +330,6 @@ namespace mu2e { } // time range of the hits auto trange = Mu2eKinKal::timeBounds(hseed.hits()); - // Begin constructing the track fit // construct the seed trajectory KTRAJ seedtraj = makeSeedTraj(hseed,trange,bnom,charge); @@ -450,23 +362,22 @@ namespace mu2e { if(!ktrk) // check that the track exists throw cet::exception("RECO")<<"mu2e::LoopHelixFit: Track fit was performed but no track is found\n"; - auto goodfit = goodFit(*ktrk); + auto goodfit = goodFit(*ktrk,seedtraj); if(print_>0) printf("[LoopHelixFit::%s] Before extending the fit: goodFit = %o, fitcon = %.4f, nHits = %2lu, %lu calo-hits\n", __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 ); - goodfit = goodFit(*ktrk); - // if finaling, apply that now. + goodfit = goodFit(*ktrk,seedtraj); + // if finalizing, apply that now. if(goodfit && fconfig_.schedule().size() > 0){ ktrk->extend(fconfig_,nohits,noexings); - goodfit = goodFit(*ktrk); + goodfit = goodFit(*ktrk,seedtraj); } } - - //store the fit quality result if it's a good fit if(print_>0) printf("[LoopHelixFit::%s] After extending the fit : goodFit = %o, fitcon = %.4f, nHits = %2lu, %lu calo-hits\n", __func__, goodfit, ktrk->fitStatus().chisq_.probability(), ktrk->strawHits().size(), ktrk->caloHits().size()); + if((!goodfit) && (! saveall_)) ktrk.reset(); return ktrk; } @@ -490,7 +401,6 @@ namespace mu2e { for(size_t iseed=0; iseed < hseedcol.size(); ++iseed) { ++nSeen_; auto const& hseed = hseedcol[iseed]; - // determine the fit direction hypotheses + check whether helix momentum is over the minimum momentum threshold (if set) auto helix_dirs = chooseHelixDir(hseed); if(helix_dirs.empty()) { @@ -498,34 +408,23 @@ namespace mu2e { continue; //bad helix, no fits to perform } ++nFit_; - const unsigned dirs_size = helix_dirs.size(); const bool undefined_dir = dirs_size > 1; //fitting multiple hypotheses to determine the best fit if(undefined_dir) ++nAmbiguous_; - // fit each track hypothesis for(auto helix_dir : helix_dirs) { auto ktrk = fitTrack(event, hseed, TrkFitDirection(helix_dir), fpart_); - if(!ktrk) continue; //ensure that the track exists - - // Check the fit - auto goodfit = goodFit(*ktrk); - // extrapolate as required - if(goodfit && extrapolate_) extrapolate(*ktrk); + if(extrap_)extrap_->extrapolate(*ktrk); if(print_>1) ktrk->printFit(std::cout,print_-1); - // save the fit result - if(goodfit || saveall_){ auto hptr = HPtr(hseedcol_h,iseed); TrkFitFlag fitflag(hptr->status()); fitflag.merge(fitflag_); - if(goodfit) fitflag.merge(TrkFitFlag::FitOK); - else fitflag.clear(TrkFitFlag::FitOK); if(undefined_dir) fitflag.merge(TrkFitFlag::AmbFitDir); // sample the fit as requested - sampleFit(*ktrk); + kkfit_.sampleFit(*ktrk); // convert to seed output format auto kkseed = kkfit_.createSeed(*ktrk,fitflag,*calo_h,*nominalTracker_h); if(print_>0) print_track_info(kkseed, *ktrk); @@ -538,7 +437,6 @@ namespace mu2e { //increment the counts if(helix_dir == TrkFitDirection::FitDirection::downstream) ++nDownstream_; if(helix_dir == TrkFitDirection::FitDirection::upstream ) ++nUpstream_; - } } //end track fit result loop } //end helix seed loop } //end helix colllection loop @@ -574,17 +472,20 @@ namespace mu2e { return ktraj; } - bool LoopHelixFit::goodFit(KKTRK const& ktrk) const { - // require physical consistency: fit can succeed but the result can have changed charge or helicity - bool retval = ktrk.fitStatus().usable() && - ktrk.fitTraj().front().parameterSign()*ktrk.seedTraj().front().parameterSign() > 0 && - ktrk.fitTraj().front().helicity()*ktrk.seedTraj().front().helicity() > 0; - // also check that the fit is inside the physical detector volume. Test where the StrawHits are + bool LoopHelixFit::goodFit(KKTRK const& ktrk,KTRAJ const& seed) const { + bool retval = ktrk.fitStatus().usable(); if(retval){ - for(auto const& shptr : ktrk.strawHits()) { - if(shptr->active() && !Mu2eKinKal::inDetector(ktrk.fitTraj().position3(shptr->time()))){ - retval = false; - break; + // require physical consistency: fit can succeed but the result can have changed charge or helicity. Test at the t0 segment + auto t0 = Mu2eKinKal::zTime(ktrk.fitTraj(),0.0,ktrk.fitTraj().range().mid()); + auto const& t0seg = ktrk.fitTraj().nearestPiece(t0); + bool retval = ktrk.fitStatus().usable() && t0seg.parameterSign()*seed.parameterSign() > 0 && t0seg.helicity()*seed.helicity() > 0; + // also check that the fit is inside the physical detector volume. Test where the StrawHits are + if(retval){ + for(auto const& shptr : ktrk.strawHits()) { + if(shptr->active() && !Mu2eKinKal::inDetector(ktrk.fitTraj().position3(shptr->time()))){ + retval = false; + break; + } } } } @@ -627,239 +528,6 @@ namespace mu2e { return true; } - void LoopHelixFit::extrapolate(KKTRK& ktrk) const { - // define the time direction according to the fit direction inside the tracker - auto const& ftraj = ktrk.fitTraj(); - if(toTrackerEnds_)toTrackerEnds(ktrk); - if(upstream_){ - auto dir0 = ftraj.direction(ftraj.t0()); - TimeDir tdir = (dir0.Z() > 0) ? TimeDir::backwards : TimeDir::forwards; - double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); - // extrapolate through the IPA in this direction. - bool exitsIPA = extrapolateIPA(ktrk,tdir); - if(exitsIPA){ // if it exits out the back, extrapolate through the target - bool exitsST = extrapolateST(ktrk,tdir); - if(exitsST) { // if it exits out the back, extrapolate to the TSDA (DS rear absorber) - bool hitTSDA = extrapolateTSDA(ktrk,tdir); - // if we hit the TSDA we are done. Otherwise if we reflected, go back through the ST - if(!hitTSDA){ // reflection upstream of the target: go back through the target - extrapolateST(ktrk,tdir); - if(backToTracker_){ // optionally extrapolate back through the IPA, then to the tracker entrance - extrapolateIPA(ktrk,tdir); - extrapolateTracker(ktrk,tdir); - } - } - } else { // reflection inside the ST; extrapolate back through the IPA, then to the tracker entrance - if(backToTracker_){ - extrapolateIPA(ktrk,tdir); - extrapolateTracker(ktrk,tdir); - } - } - } else { // reflection inside the IPA; extrapolate back through the IPA, then to the tracker entrance - if(backToTracker_)ktrk.extrapolate(tdir,trackerFront_); - } - // optionally test for intersection with the OPA - if(toOPA_)toOPA(ktrk,starttime,tdir); - } - } - - void LoopHelixFit::toTrackerEnds(KKTRK& ktrk) const { - // 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 volumei - auto const& ftraj = ktrk.fitTraj(); - auto dir0 = ftraj.direction(ftraj.t0()); - TimeDir fronttdir = (dir0.Z() > 0) ? TimeDir::backwards : TimeDir::forwards; - TimeDir backtdir = (dir0.Z() > 0) ? TimeDir::forwards : TimeDir::backwards; - auto tofront = ktrk.extrapolate(fronttdir,trackerFront_); - auto toback = ktrk.extrapolate(backtdir,trackerBack_); - // record the standard tracker intersections - static const SurfaceId tt_front("TT_Front"); - static const SurfaceId tt_mid("TT_Mid"); - static const SurfaceId tt_back("TT_Back"); - - // start with the middle - auto midinter = KinKal::intersect(ftraj,*trkmidptr_,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,*trkfrontptr_,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,*trkfrontptr_,frange,intertol_,fronttdir); - } - if(frontinter.good()) ktrk.addIntersection(tt_front,frontinter); - } - if(toback){ - // 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,*trkbackptr_,brange,intertol_,backtdir); - if(backinter.good())ktrk.addIntersection(tt_back,backinter); - } - } - - bool LoopHelixFit::extrapolateIPA(KKTRK& ktrk,TimeDir tdir) const { - if(extrapIPA_.debug() > 0)std::cout << "extrapolating to IPA " << std::endl; - // extraplate the fit through the IPA. This will add material effects for each intersection. It will continue till the - // track exits the IPA - extrapIPA_.reset(); - auto const& ftraj = ktrk.fitTraj(); - static const SurfaceId IPASID("IPA"); - double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); - auto startdir = ftraj.direction(starttime); - do { - ktrk.extrapolate(tdir,extrapIPA_); - 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 = smap_.DS().innerProtonAbsorberPtr(); - KKIPAXINGPTR ipaxingptr = std::make_shared(IPA,IPASID,*kkmat_.IPAMaterial(),extrapIPA_.intersection(),reftrajptr,ipathick_,extrapIPA_.interTolerance()); - if(extrapIPA_.debug() > 0){ - double dmom, paramomvar, perpmomvar; - ipaxingptr->materialEffects(dmom,paramomvar,perpmomvar); - std::cout << "IPA Xing dmom " << dmom << " para momsig " << sqrt(paramomvar) << " perp momsig " << sqrt(perpmomvar) << std::endl; - std::cout << " before append mom = " << reftrajptr->momentum(); - } - ktrk.addIPAXing(ipaxingptr,tdir); - if(extrapIPA_.debug() > 0){ - auto const& newtrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); - std::cout << " after append mom = " << newtrajptr->momentum() << std::endl; - } - } - } while(extrapIPA_.intersection().good()); - // check if the particle exited in the same physical direction or not (reflection) - double endtime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); - auto enddir = ftraj.direction(endtime); - if(enddir.Z() * startdir.Z() > 0.0){ - return true; - } - return false; - } - - bool LoopHelixFit::extrapolateST(KKTRK& ktrk,TimeDir tdir) const { - // 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 - extrapST_.reset(); - auto const& ftraj = ktrk.fitTraj(); - double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); - auto startdir = ftraj.direction(starttime); - if(extrapST_.debug() > 0)std::cout << "extrapolating to ST " << std::endl; - do { - ktrk.extrapolate(tdir,extrapST_); - 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()); - if(extrapST_.debug() > 0){ - double dmom, paramomvar, perpmomvar; - stxingptr->materialEffects(dmom,paramomvar,perpmomvar); - std::cout << "ST Xing dmom " << dmom << " para momsig " << sqrt(paramomvar) << " perp momsig " << sqrt(perpmomvar) << std::endl; - std::cout << " before append mom = " << reftrajptr->momentum(); - } - // Add the xing. This truncates the fit - ktrk.addSTXing(stxingptr,tdir); - if(extrapST_.debug() > 0){ - auto const& newtrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); - std::cout << " after append mom = " << newtrajptr->momentum() << std::endl; - } - } - } while(extrapST_.intersection().good()); - // check if the particle exited in the same physical direction or not (reflection) - double endtime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); - auto enddir = ftraj.direction(endtime); - if(enddir.Z() * startdir.Z() > 0.0){ - return true; - } - return false; - } - - bool LoopHelixFit::extrapolateTracker(KKTRK& ktrk,TimeDir tdir) const { - 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,*trkfrontptr_,ktraj.range(),intertol_,tdir); - if(trkfrontinter.onsurface_){ // dont worry about bounds here - ktrk.addIntersection(TrackerSID,trkfrontinter); - return true; - } - return false; - } - - bool LoopHelixFit::extrapolateTSDA(KKTRK& ktrk,TimeDir tdir) const { - if(TSDA_.debug() > 0)std::cout << "extrapolating to TSDA " << std::endl; - auto const& ftraj = ktrk.fitTraj(); - static const SurfaceId TSDASID("TSDA"); - ktrk.extrapolate(tdir,TSDA_); - // if we reflected we're done. Otherwize, save the TSDA intersection - double tend = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); - auto epos = ftraj.position3(tend); - bool retval = epos.Z() < TSDA_.zVal(); - if(retval){ - auto const& ktraj = tdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); - auto tsdainter = KinKal::intersect(ftraj,*tsdaptr_,ktraj.range(),intertol_,tdir); - if(tsdainter.onsurface_)ktrk.addIntersection(TSDASID,tsdainter); - } - return retval; - } - - void LoopHelixFit::toOPA(KKTRK& ktrk, double tstart, TimeDir tdir) const { - 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,*opaptr_,trange,intertol_,tdir); - if(opainter.good()){ - ktrk.addIntersection(OPASID,opainter); - } - } - - void LoopHelixFit::sampleFit(KKTRK& kktrk) const { - auto const& ftraj = kktrk.fitTraj(); - std::vector ranges; - // test for reflection, and if present, split the test in 2 - auto refltraj = ftraj.reflection(ftraj.range().begin()); - if(refltraj){ - double tmid = refltraj->range().begin(); - ranges.push_back(TimeRange(ftraj.range().begin(),tmid)); - ranges.push_back(TimeRange(tmid,ftraj.range().end())); - } else { - ranges.push_back(ftraj.range()); - } - for(auto range : ranges) { - double tbeg = range.begin(); - double tend = range.end(); - for(auto const& surf : sample_){ - // search for intersections with each surface within the specified time range, going forwards in time - bool goodinter(true); - size_t max_inter = 100; // limit the number of intersections - size_t cur_inter = 0; - // loop to find multiple intersections - while(goodinter && tbeg < tend && cur_inter < max_inter){ - cur_inter += 1; - TimeRange irange(tbeg,tend); - auto surfinter = KinKal::intersect(ftraj,*surf.second,irange,intertol_); - goodinter = surfinter.onsurface_ && ( (! sampleinbounds_) || surfinter.inbounds_ ) && ( (!sampleinrange_) || surfinter.inrange_); - if(goodinter) { - // save the intersection information - kktrk.addIntersection(surf.first,surfinter); - // update for the next intersection - // move past existing intersection to avoid repeating - double epsilon = intertol_/ftraj.speed(surfinter.time_); - tbeg = surfinter.time_ + epsilon; - } - if(print_>1) printf(" [LoopHelixFit::%s::%s] Found intersection with surface %15s\n", - __func__, moduleDescription().moduleLabel().c_str(), surf.first.name().c_str()); - } - } - } - } - void LoopHelixFit::print_track_info(const KalSeed& kkseed, const KKTRK& ktrk) const { if(print_>0) printf("[LoopHelixFit::%s::%s] Create seed : fitcon = %.4f, nHits = %2lu, seedActiveHits = %2u, %lu calo-hits\n", __func__, moduleDescription().moduleLabel().c_str(), ktrk.fitStatus().chisq_.probability(), ktrk.strawHits().size(), diff --git a/Mu2eKinKal/src/RegrowKinematicLine_module.cc b/Mu2eKinKal/src/RegrowKinematicLine_module.cc new file mode 100644 index 0000000000..8b872cb1cc --- /dev/null +++ b/Mu2eKinKal/src/RegrowKinematicLine_module.cc @@ -0,0 +1,278 @@ +// +// Instantiation of RegrowKalSeed for KinematicLine fits +// +// Original author: D. Brown (LBNL) 4/18/2025 +// +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/Sequence.h" +#include "fhiclcpp/types/Table.h" +#include "fhiclcpp/types/OptionalTable.h" +#include "fhiclcpp/types/Tuple.h" +#include "fhiclcpp/types/OptionalAtom.h" +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/Handle.h" +// conditions +#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh" +#include "Offline/ProditionsService/inc/ProditionsHandle.hh" +#include "Offline/TrackerConditions/inc/StrawResponse.hh" +#include "Offline/BFieldGeom/inc/BFieldManager.hh" +#include "Offline/GlobalConstantsService/inc/ParticleDataList.hh" +#include "Offline/DataProducts/inc/SurfaceId.hh" +#include "Offline/KinKalGeom/inc/SurfaceMap.hh" +// utiliites +#include "Offline/GeometryService/inc/GeomHandle.hh" +#include "Offline/TrackerGeom/inc/Tracker.hh" +#include "Offline/GeometryService/inc/GeometryService.hh" +#include "Offline/GeneralUtilities/inc/Angles.hh" +#include "Offline/TrkReco/inc/TrkUtilities.hh" +#include "Offline/CalorimeterGeom/inc/Calorimeter.hh" +#include "Offline/GeneralUtilities/inc/OwningPointerCollection.hh" +// data +#include "Offline/DataProducts/inc/PDGCode.hh" +#include "Offline/DataProducts/inc/Helicity.hh" +#include "Offline/RecoDataProducts/inc/ComboHit.hh" +#include "Offline/RecoDataProducts/inc/StrawHitFlag.hh" +#include "Offline/RecoDataProducts/inc/KalSeed.hh" +#include "Offline/RecoDataProducts/inc/TrkFitDirection.hh" +#include "Offline/DataProducts/inc/IndexMap.hh" +#include "Offline/MCDataProducts/inc/KalSeedMC.hh" +// KinKal +#include "KinKal/Fit/Track.hh" +#include "KinKal/Fit/Config.hh" +#include "KinKal/General/Parameters.hh" +#include "KinKal/General/Vectors.hh" +#include "KinKal/Geometry/Cylinder.hh" +#include "KinKal/Geometry/Disk.hh" +#include "KinKal/Geometry/Frustrum.hh" +#include "KinKal/Trajectory/KinematicLine.hh" +#include "KinKal/Trajectory/ParticleTrajectory.hh" +#include "KinKal/Trajectory/PiecewiseClosestApproach.hh" +#include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" +// Mu2eKinKal +#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/Mu2eKinKal/inc/KKStrawHit.hh" +#include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" +#include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" +#include "Offline/Mu2eKinKal/inc/KKCaloHit.hh" +#include "Offline/Mu2eKinKal/inc/KKBField.hh" +#include "Offline/Mu2eKinKal/inc/KKExtrap.hh" +// C++ +#include +#include +#include +#include +#include + +namespace mu2e { + using KinKal::VEC3; + using KinKal::DMAT; + using KinKal::DVEC; + using KinKal::TimeDir; + using KinKal::Config; + using MatEnv::DetMaterial; + using KKConfig = Mu2eKinKal::KinKalConfig; + using Mu2eKinKal::KKFinalConfig; + using KKFitConfig = Mu2eKinKal::KKFitConfig; + using KKModuleConfig = Mu2eKinKal::KKModuleConfig; + using KKMaterialConfig = KKMaterial::Config; + using SDIS = std::set; + + using Name = fhicl::Name; + using Comment = fhicl::Comment; + struct RegrowKinematicLineConfig { + fhicl::Atom debug{Name("debug"), Comment("Debug level"), 0}; + fhicl::Atom kalSeedCollection {Name("KalSeedCollection"), Comment("KalSeed collection to process") }; + fhicl::Atom comboHitCollection {Name("ComboHitCollection"), Comment("Reduced ComboHit collection") }; + 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") }; + }; + + class RegrowKinematicLine : public art::EDProducer { + public: + using Parameters = art::EDProducer::Table; + using KTRAJ = KinKal::KinematicLine; + using PKTRAJ = KinKal::ParticleTrajectory; + using PKTRAJPTR = std::unique_ptr; + using KKTRK = KKTrack; + using KKTRKCOL = OwningPointerCollection; + using KKSTRAWHIT = KKStrawHit; + using KKSTRAWHITPTR = std::shared_ptr; + using KKSTRAWHITCOL = std::vector; + using KKSTRAWXING = KKStrawXing; + using KKSTRAWXINGPTR = std::shared_ptr; + using KKSTRAWXINGCOL = std::vector; + using KKIPAXING = KKShellXing; + using KKIPAXINGPTR = std::shared_ptr; + using KKIPAXINGCOL = std::vector; + using KKSTXING = KKShellXing; + using KKSTXINGPTR = std::shared_ptr; + using KKSTXINGCOL = std::vector; + using KKCALOHIT = KKCaloHit; + using KKCALOHITPTR = std::shared_ptr; + using KKCALOHITCOL = std::vector; + using KKFIT = KKFit; + + using MEAS = KinKal::Hit; + using MEASPTR = std::shared_ptr; + using MEASCOL = std::vector; + using EXING = KinKal::ElementXing; + using EXINGPTR = std::shared_ptr; + using EXINGCOL = std::vector; + 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; + void endJob() override; + private: + int debug_; + ProditionsHandle strawResponse_h_; + ProditionsHandle alignedTracker_h_; + 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_; + art::InputTag ksmca_T_; + bool fillMCAssns_; + std::unique_ptr extrap_; + }; + + RegrowKinematicLine::RegrowKinematicLine(const Parameters& settings) : art::EDProducer(settings), + 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())), + fillMCAssns_(settings().kalSeedMCAssns(ksmca_T_)) + { + produces(); + produces(); + if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings(),kkmat_); + if( fillMCAssns_){ + consumes(ksmca_T_); + produces (); + } + } + + void RegrowKinematicLine::beginRun(art::Run& run) + { + GeomHandle bfmgr; + GeomHandle det; + kkbf_ = std::move(std::make_unique(*bfmgr,*det)); + } + + void RegrowKinematicLine::produce(art::Event& event) + { + // proditions + auto const& strawresponse = strawResponse_h_.getPtr(event.id()); + auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); + GeomHandle nominalTracker_h; + GeomHandle calo_h; + // find input event data + auto kseed_H = event.getValidHandle(kseedcol_T_); + const auto& kseedcol = *kseed_H; + auto ch_H = event.getValidHandle(chcol_T_); + const auto& chcol = *ch_H; + auto indexmap_H = event.getValidHandle(indexmap_T_); + const auto& indexmap = *indexmap_H; + auto KalSeedCollectionPID = event.getProductID(); + auto KalSeedCollectionGetter = event.productGetter(KalSeedCollectionPID); + art::Handle ksmca_H; + // create outputs + unique_ptr ktrkcol(new KKTRKCOL ); + unique_ptr rgkseedcol(new KalSeedCollection ); + std::unique_ptr ksmca; + // deal with MC + if(fillMCAssns_){ + ksmca_H = event.getHandle(ksmca_T_); + if(!ksmca_H)throw cet::exception("RECO")<<"mu2e::RegrowKinematicLine: No KalSeedMCAssns found" << endl; + ksmca = std::unique_ptr(new KalSeedMCAssns); + } + size_t iseed(0); + for (auto const& kseed : kseedcol) { + if(!kseed.kinematicLineFit())throw cet::exception("RECO")<<"mu2e::RegrowKinematicLine: passed KalSeed from non-KinematicLine fit " << endl; + // regrow the components from the seed + PKTRAJPTR trajptr = kseed.kinematicLineFitTrajectory(); + KKSTRAWHITCOL strawhits; + strawhits.reserve(kseed.hits().size()); + KKSTRAWXINGCOL strawxings; + strawxings.reserve(kseed.straws().size()); + KKCALOHITCOL calohits; + 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(), + 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; + } + if(debug_ > 2){ + unsigned nhactive(0); + unsigned nsactive(0); + for( auto const& strawh : strawhits)if(strawh->active())++nhactive; + for( auto const& strawx : strawxings)if(strawx->active())++nsactive; + std::cout << "Regrew " << nhactive << " active hits and " << nsactive << " active straws" << std::endl; + + } + if(debug_ > 5)static_cast*>(trajptr.get())->print(std::cout,2); + if(goodhits){ + // create the KKTrack from these components + auto ktrk = std::make_unique(config_,*kkbf_,kseed.particle(),trajptr,strawhits,strawxings,calohits,domains); + if(ktrk && ktrk->fitStatus().usable()){ + if(debug_ > 0) std::cout << "RegrowKinematicLine: successful track refit" << std::endl; + // extrapolate as requested + if(extrap_)extrap_->extrapolate(*ktrk); + // sample the fit as requested + kkfit_.sampleFit(*ktrk); + // convert to seed output format + TrkFitFlag fitflag = kseed.status(); + fitflag.merge(TrkFitFlag::Regrown); + auto rgks = kkfit_.createSeed(*ktrk,fitflag,*calo_h,*nominalTracker_h); + rgkseedcol->push_back(rgks); + if(fillMCAssns_){ + // find the MC assns + auto ksmcai = (*ksmca_H)[iseed]; + auto origksp = art::Ptr(kseed_H,iseed); + // test this is the right ptr + if(ksmcai.first != origksp)throw cet::exception("Reco")<<"mu2e::RegrowKinematicLine: wrong KalSeed ptr"<< std::endl; + auto mcseedp = ksmcai.second; + auto rgksp = art::Ptr(KalSeedCollectionPID,rgkseedcol->size()-1,KalSeedCollectionGetter); + ksmca->addSingle(rgksp,mcseedp); + // add the original too + ksmca->addSingle(origksp,mcseedp); + } + if(debug_ > 5)static_cast&>(ktrk->fitTraj()).print(std::cout,2); + ktrkcol->push_back(ktrk.release()); + } else if(debug_ > 0) + std::cout << "RegrowKinematicLine: failed track refit, status " << ktrk->fitStatus() << std::endl; + } + ++iseed; + } + // store output + event.put(move(ktrkcol)); + event.put(move(rgkseedcol)); + if(fillMCAssns_)event.put(move(ksmca)); + } + + void RegrowKinematicLine::endJob() + {} + +} +DEFINE_ART_MODULE(mu2e::RegrowKinematicLine) diff --git a/Mu2eKinKal/src/RegrowLoopHelix_module.cc b/Mu2eKinKal/src/RegrowLoopHelix_module.cc new file mode 100644 index 0000000000..ffd1d884ff --- /dev/null +++ b/Mu2eKinKal/src/RegrowLoopHelix_module.cc @@ -0,0 +1,289 @@ +// +// Instantiation of RegrowKalSeed for LoopHelix fits +// +// Original author: D. Brown (LBNL) 4/18/2025 +// +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/Sequence.h" +#include "fhiclcpp/types/Table.h" +#include "fhiclcpp/types/OptionalTable.h" +#include "fhiclcpp/types/Tuple.h" +#include "fhiclcpp/types/OptionalAtom.h" +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/Handle.h" +// conditions +#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh" +#include "Offline/ProditionsService/inc/ProditionsHandle.hh" +#include "Offline/TrackerConditions/inc/StrawResponse.hh" +#include "Offline/BFieldGeom/inc/BFieldManager.hh" +#include "Offline/GlobalConstantsService/inc/ParticleDataList.hh" +#include "Offline/DataProducts/inc/SurfaceId.hh" +#include "Offline/KinKalGeom/inc/SurfaceMap.hh" +// utiliites +#include "Offline/GeometryService/inc/GeomHandle.hh" +#include "Offline/TrackerGeom/inc/Tracker.hh" +#include "Offline/GeometryService/inc/GeometryService.hh" +#include "Offline/GeneralUtilities/inc/Angles.hh" +#include "Offline/TrkReco/inc/TrkUtilities.hh" +#include "Offline/CalorimeterGeom/inc/Calorimeter.hh" +#include "Offline/GeneralUtilities/inc/OwningPointerCollection.hh" +// data +#include "Offline/DataProducts/inc/PDGCode.hh" +#include "Offline/DataProducts/inc/Helicity.hh" +#include "Offline/RecoDataProducts/inc/ComboHit.hh" +#include "Offline/RecoDataProducts/inc/StrawHitFlag.hh" +#include "Offline/RecoDataProducts/inc/KalSeed.hh" +#include "Offline/RecoDataProducts/inc/TrkFitDirection.hh" +#include "Offline/DataProducts/inc/IndexMap.hh" +#include "Offline/MCDataProducts/inc/KalSeedMC.hh" +// KinKal +#include "KinKal/Fit/Track.hh" +#include "KinKal/Fit/Config.hh" +#include "KinKal/General/Parameters.hh" +#include "KinKal/General/Vectors.hh" +#include "KinKal/Geometry/Cylinder.hh" +#include "KinKal/Geometry/Disk.hh" +#include "KinKal/Geometry/Frustrum.hh" +#include "KinKal/Trajectory/LoopHelix.hh" +#include "KinKal/Trajectory/ParticleTrajectory.hh" +#include "KinKal/Trajectory/PiecewiseClosestApproach.hh" +#include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" +// Mu2eKinKal +#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/Mu2eKinKal/inc/KKStrawHit.hh" +#include "Offline/Mu2eKinKal/inc/KKStrawHitCluster.hh" +#include "Offline/Mu2eKinKal/inc/KKStrawXing.hh" +#include "Offline/Mu2eKinKal/inc/KKCaloHit.hh" +#include "Offline/Mu2eKinKal/inc/KKBField.hh" +#include "Offline/Mu2eKinKal/inc/KKExtrap.hh" +// C++ +#include +#include +#include +#include +#include + +namespace mu2e { + using KinKal::VEC3; + using KinKal::DMAT; + using KinKal::DVEC; + using KinKal::TimeDir; + using KinKal::Config; + using MatEnv::DetMaterial; + using KKConfig = Mu2eKinKal::KinKalConfig; + using Mu2eKinKal::KKFinalConfig; + using KKFitConfig = Mu2eKinKal::KKFitConfig; + using KKModuleConfig = Mu2eKinKal::KKModuleConfig; + using KKMaterialConfig = KKMaterial::Config; + using SDIS = std::set; + + using Name = fhicl::Name; + using Comment = fhicl::Comment; + struct RegrowLoopHelixConfig { + fhicl::Atom debug{Name("debug"), Comment("Debug level"), 0}; + fhicl::Atom kalSeedCollection {Name("KalSeedCollection"), Comment("KalSeed collection to process") }; + fhicl::Atom comboHitCollection {Name("ComboHitCollection"), Comment("Reduced ComboHit collection") }; + fhicl::Atom caloClusterCollection {Name("CaloClusterCollection"), Comment("CaloCluster collection ") }; + 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") }; + + fhicl::OptionalTable extrapSettings { Name("ExtrapolationSettings") }; + }; + + class RegrowLoopHelix : public art::EDProducer { + public: + using Parameters = art::EDProducer::Table; + using KTRAJ = KinKal::LoopHelix; + using PKTRAJ = KinKal::ParticleTrajectory; + using PKTRAJPTR = std::unique_ptr; + using KKTRK = KKTrack; + using KKTRKCOL = OwningPointerCollection; + using KKSTRAWHIT = KKStrawHit; + using KKSTRAWHITPTR = std::shared_ptr; + using KKSTRAWHITCOL = std::vector; + using KKSTRAWXING = KKStrawXing; + using KKSTRAWXINGPTR = std::shared_ptr; + using KKSTRAWXINGCOL = std::vector; + using KKIPAXING = KKShellXing; + using KKIPAXINGPTR = std::shared_ptr; + using KKIPAXINGCOL = std::vector; + using KKSTXING = KKShellXing; + using KKSTXINGPTR = std::shared_ptr; + using KKSTXINGCOL = std::vector; + using KKCALOHIT = KKCaloHit; + using KKCALOHITPTR = std::shared_ptr; + using KKCALOHITCOL = std::vector; + using KKFIT = KKFit; + + using MEAS = KinKal::Hit; + using MEASPTR = std::shared_ptr; + using MEASCOL = std::vector; + using EXING = KinKal::ElementXing; + using EXINGPTR = std::shared_ptr; + using EXINGCOL = std::vector; + 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; + void endJob() override; + private: + int debug_; + ProditionsHandle strawResponse_h_; + ProditionsHandle alignedTracker_h_; + 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_; + art::ProductToken indexmap_T_; + art::InputTag ksmca_T_; + bool fillMCAssns_; + bool extend_; + std::unique_ptr extrap_; + }; + + RegrowLoopHelix::RegrowLoopHelix(const Parameters& settings) : art::EDProducer(settings), + 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())), + indexmap_T_(consumes(settings().indexMap())), + fillMCAssns_(settings().kalSeedMCAssns(ksmca_T_)), + extend_(settings().extend()) + { + produces(); + produces(); + if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings(),kkmat_); + if( fillMCAssns_){ + consumes(ksmca_T_); + produces (); + } + } + + void RegrowLoopHelix::beginRun(art::Run& run) + { + GeomHandle bfmgr; + GeomHandle det; + kkbf_ = std::move(std::make_unique(*bfmgr,*det)); + } + + void RegrowLoopHelix::produce(art::Event& event) + { + // proditions + auto const& strawresponse = strawResponse_h_.getPtr(event.id()); + auto const& tracker = alignedTracker_h_.getPtr(event.id()).get(); + GeomHandle nominalTracker_h; + GeomHandle calo_h; + // find input event data + auto kseed_H = event.getValidHandle(kseedcol_T_); + const auto& kseedcol = *kseed_H; + auto ch_H = event.getValidHandle(chcol_T_); + const auto& chcol = *ch_H; + auto cc_H = event.getHandle(cccol_T_); + auto indexmap_H = event.getValidHandle(indexmap_T_); + const auto& indexmap = *indexmap_H; + auto KalSeedCollectionPID = event.getProductID(); + auto KalSeedCollectionGetter = event.productGetter(KalSeedCollectionPID); + art::Handle ksmca_H; + // create outputs + unique_ptr ktrkcol(new KKTRKCOL ); + unique_ptr rgkseedcol(new KalSeedCollection ); + std::unique_ptr ksmca; + // deal with MC + if(fillMCAssns_){ + ksmca_H = event.getHandle(ksmca_T_); + if(!ksmca_H)throw cet::exception("RECO")<<"mu2e::RegrowLoopHelix: No KalSeedMCAssns found" << endl; + ksmca = std::unique_ptr(new KalSeedMCAssns); + } + size_t iseed(0); + for (auto const& kseed : kseedcol) { + if(!kseed.loopHelixFit())throw cet::exception("RECO")<<"mu2e::RegrowLoopHelix: passed KalSeed from non-LoopHelix fit " << endl; + // regrow the components from the seed + PKTRAJPTR trajptr = kseed.loopHelixFitTrajectory(); + KKSTRAWHITCOL strawhits; + strawhits.reserve(kseed.hits().size()); + KKSTRAWXINGCOL strawxings; + strawxings.reserve(kseed.straws().size()); + KKCALOHITCOL calohits; + 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(), + 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; + } + if(debug_ > 2){ + unsigned nhactive(0); + unsigned nsactive(0); + for( auto const& strawh : strawhits)if(strawh->active())++nhactive; + for( auto const& strawx : strawxings)if(strawx->active())++nsactive; + std::cout << "Regrew " << nhactive << " active hits and " << nsactive << " active straws" << std::endl; + } + if(debug_ > 5)static_cast*>(trajptr.get())->print(std::cout,2); + // require hits and consistent BField domains + if(goodhits && (domains.size() > 0 || !config_.bfcorr_)){ + // create the KKTrack from these components + 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(ktrk->fitStatus().usable()){ + // extrapolate as requested + if(extrap_)extrap_->extrapolate(*ktrk); + // sample the fit as requested + kkfit_.sampleFit(*ktrk); + // convert to seed output format + TrkFitFlag fitflag = kseed.status(); + fitflag.merge(TrkFitFlag::Regrown); + auto rgks = kkfit_.createSeed(*ktrk,fitflag,*calo_h,*nominalTracker_h); + rgkseedcol->push_back(rgks); + if(fillMCAssns_){ + // find the MC assns + auto ksmcai = (*ksmca_H)[iseed]; + auto origksp = art::Ptr(kseed_H,iseed); + // test this is the right ptr + if(ksmcai.first != origksp)throw cet::exception("Reco")<<"mu2e::RegrowLoopHelix: wrong KalSeed ptr"<< std::endl; + auto mcseedp = ksmcai.second; + auto rgksp = art::Ptr(KalSeedCollectionPID,rgkseedcol->size()-1,KalSeedCollectionGetter); + ksmca->addSingle(rgksp,mcseedp); + // add the original too + ksmca->addSingle(origksp,mcseedp); + } + if(debug_ > 5)static_cast&>(ktrk->fitTraj()).print(std::cout,2); + ktrkcol->push_back(ktrk.release()); + } + } else if(debug_ > 0) + std::cout << "RegrowLoopHelix: failed track refit, status " << ktrk->fitStatus() << std::endl; + } + ++iseed; + } + // store output + event.put(move(ktrkcol)); + event.put(move(rgkseedcol)); + if(fillMCAssns_)event.put(move(ksmca)); + } + + void RegrowLoopHelix::endJob() + {} + +} +DEFINE_ART_MODULE(mu2e::RegrowLoopHelix) diff --git a/Print/src/KalSeedPrinter.cc b/Print/src/KalSeedPrinter.cc index f15517c2e8..3b76f7b4c3 100644 --- a/Print/src/KalSeedPrinter.cc +++ b/Print/src/KalSeedPrinter.cc @@ -70,11 +70,13 @@ void mu2e::KalSeedPrinter::Print(const mu2e::KalSeed& obj, int ind, if(obj.loopHelixFit())fittype = "LoopHelix"; if(obj.centralHelixFit())fittype = "CentralHelix"; if(obj.kinematicLineFit())fittype = "KinematicLine"; - os << std::setprecision(2) << obj.particle() << " " << std::setw(2) + unsigned nactive(0); + for(auto const& hit : obj.hits())if(hit._ambig>WireHitState::inactive)nactive++; + os << " " << std::setprecision(2) << obj.particle() << " " << std::setw(2) << std::setprecision(2) << fittype << " " << std::setw(2) << std::setprecision(3) << obj.fitConsistency() << " " << std::setw(5) << std::setprecision(1) << (obj.caloCluster().isNull() ? "no" : "yes") << std::setw(3) - << std::setprecision(4) << obj.hits().size() << " " << std::setw(3) + << std::setprecision(4) << nactive << " " << std::setw(3) << std::setprecision(4) << obj.straws().size() << " " << std::setw(3) << std::setprecision(4) << obj.segments().size() << " " << std::setw(3) << std::setprecision(4) << obj.intersections().size() << " " << std::setw(3); @@ -87,10 +89,17 @@ void mu2e::KalSeedPrinter::Print(const mu2e::KalSeed& obj, int ind, } os << "\n"; - } else if (verbose() >= 2) { - os << " intersections: \n"; + } + if(verbose() >= 2) { + os << std::setprecision(3) << "Intersections: \n"; for (auto const& inter : obj.intersections()) { - os << " sid " << inter.surfaceId() << " time " << inter.time() << " pos "<< inter.position3() << " P " << inter.momentum3() << " dP " << inter.dMom() << "\n"; + os << " sid " << inter.surfaceId() << " time " << inter.time() << " Pos " << inter.position3() << " Mom " << inter.momentum3() << " dP " << inter.dMom() << "\n"; + } + } + if(verbose() >= 3) { + os << std::setprecision(4) << "Segments: \n"; + for (auto const& iseg : obj.segments()) { + os << "Range" << iseg.timeRange() << " Pos " << iseg.position3() << " Mom " << iseg.momentum3() << " Bnom " << iseg.bnom() << "\n"; } } } diff --git a/RecoDataProducts/inc/KalIntersection.hh b/RecoDataProducts/inc/KalIntersection.hh index 92b88eddf1..0b2ec877fa 100644 --- a/RecoDataProducts/inc/KalIntersection.hh +++ b/RecoDataProducts/inc/KalIntersection.hh @@ -14,7 +14,7 @@ #include namespace mu2e { struct KalIntersection { - KinKal::ParticleStateEstimate pstate_; // particle state at intersection point/time + KinKal::ParticleStateEstimate pstate_; // particle state at intersection point/time before traversing any material XYZVectorF bnom_; // Bfield at this intersection, needed to reconstitute trajectory SurfaceId surfid_; // which surface in the reco geometry was interestected KinKal::Intersection kkinter_; // kinkal intersection diff --git a/RecoDataProducts/inc/KalSeed.hh b/RecoDataProducts/inc/KalSeed.hh index e7e71a093e..56f72e9ed9 100644 --- a/RecoDataProducts/inc/KalSeed.hh +++ b/RecoDataProducts/inc/KalSeed.hh @@ -48,10 +48,12 @@ namespace mu2e { auto const& segments() const { return _segments; } auto const& intersections() const { return _inters; } auto const& status() const { return _status; } + auto const& domainBounds() const { return _domainbounds; } double t0Val() const; double t0Var() const; float chisquared() const { return _chisq; } int nDOF() const { return _ndof; } + float domainTolerance() const { return _dtol; } unsigned nHits(bool active=true) const; float fitConsistency() const { return _fitcon; } UInt_t nTrajSegments() const { return _segments.size(); } @@ -72,22 +74,26 @@ namespace mu2e { KLPTPtr kinematicLineFitTrajectory() const; // global information about the track - PDGCode::type _tpart = PDGCode::unknown; // particle assumed for this fit - TrkFitFlag _status; // status of this fit: includes alglorithm information + PDGCode::type _tpart = PDGCode::unknown; // particle assumed for this fit + TrkFitFlag _status; // status of this fit: includes alglorithm information float _chisq = -1; // fit chisquared value - int _ndof = -1; + int _ndof = -1; float _fitcon = -1; // fit consistency float _maxgap = 0; float _avggap = 0; // information about trajectory gaps + float _dtol = 0; // tolerance used when creating BField domain // // contained content substructure. // - std::vector _segments; // segments of the Kalman filter fit result - std::vector _inters; // intersections with materials or reference locations - std::vector _hits; // hit seeds for all the hits used in this fit - std::vector _hitcalibs; // extra calibration/alignment info - std::vector _straws; // straws interesected by this fit - TrkCaloHitSeed _chit; // CaloCluster-based hit. If it has no CaloCluster, this has no content + std::vector _segments; // segments of the Kalman filter fit result + std::vector _inters; // intersections with materials or reference locations + std::vector _hits; // hit seeds for all the hits used in this fit + std::vector _hitcalibs; // extra calibration/alignment info + std::vector _straws; // straws interesected by this fit + std::vector _domainbounds; // domain time boundaries + TrkCaloHitSeed _chit; // CaloCluster-based hit. If it has no CaloCluster, this has no content + // static value used in regrowing + static const double _regrowtol; // Minimimum time length for adding a segment // // deprecated BTrk legacy content, DO NOT write any new code which depends on these functions // find the nearest segment to a given the time @@ -95,7 +101,7 @@ namespace mu2e { std::vector::const_iterator nearestSegment(const XYZVectorF& pos) const; // find nearest segment to a GLOBAL position float flt0() const { return _flt0; } HitT0 t0() const; - float _flt0 = 0.0; // flight distance where the track crosses the tracker midplane (z=0). Redundant with t0 in KinKal fits, and in the wrong unit + float _flt0 = 0.0; // flight distance where the track crosses the tracker midplane (z=0). Redundant with t0 in KinKal fits, and in the wrong unit }; typedef std::vector KalSeedCollection; typedef art::Ptr KalSeedPtr; diff --git a/RecoDataProducts/inc/KalSeedAssns.hh b/RecoDataProducts/inc/KalSeedAssns.hh index e34fb8ab94..467b70f45c 100644 --- a/RecoDataProducts/inc/KalSeedAssns.hh +++ b/RecoDataProducts/inc/KalSeedAssns.hh @@ -5,7 +5,7 @@ #include "canvas/Persistency/Common/Assns.h" namespace mu2e { - // Assns between a KalSeed and the HelixSeed it was created from + // Assns between a KalSeed and the thing it was created from typedef art::Assns KalHelixAssns; typedef art::Assns KalLineAssns; } diff --git a/RecoDataProducts/inc/TrkFitFlag.hh b/RecoDataProducts/inc/TrkFitFlag.hh index b49442680a..e0d7c8126a 100644 --- a/RecoDataProducts/inc/TrkFitFlag.hh +++ b/RecoDataProducts/inc/TrkFitFlag.hh @@ -19,6 +19,7 @@ namespace mu2e { circleConverged,phizConverged,helixConverged,seedConverged,kalmanConverged, //12 MatCorr, BFCorr, FitOK, //15 KSF=16, KFF, TPRHelix, CPRHelix, Straight, KKLoopHelix,KKCentralHelix,KKLine, APRHelix, MPRHelix, AmbFitDir, //26 + Regrown=27, MCSeed=31, //31 }; diff --git a/RecoDataProducts/inc/TrkStraw.hh b/RecoDataProducts/inc/TrkStraw.hh index 8f48e024cd..3af028030b 100644 --- a/RecoDataProducts/inc/TrkStraw.hh +++ b/RecoDataProducts/inc/TrkStraw.hh @@ -19,6 +19,7 @@ namespace mu2e { _poca(pocadata.sensorPoca().Vect()), _doca(pocadata.doca()), _docavar(pocadata.docaVar()), + _toca(pocadata.particleToca()), _dirdot(pocadata.dirDot()), _radlen(radlen), _dmom(dmom) @@ -41,6 +42,7 @@ namespace mu2e { XYZVectorF _poca; // POCA to the straw axis float _doca = 0.0; // DOCA from the track to the straw axis float _docavar = 0.0; // DOCA variance + float _toca = 0.0; // TOCA from the track to the straw axis float _dirdot = 0.0; // dot product between straw axis and track direction float _gaspath = 0.0; // path length in gas material float _wallpath = 0.0; // path length in straw wall material diff --git a/RecoDataProducts/inc/TrkStrawHitSeed.hh b/RecoDataProducts/inc/TrkStrawHitSeed.hh index 2afbfa407f..dc56783c9c 100644 --- a/RecoDataProducts/inc/TrkStrawHitSeed.hh +++ b/RecoDataProducts/inc/TrkStrawHitSeed.hh @@ -8,6 +8,7 @@ #include "Offline/RecoDataProducts/inc/HitT0.hh" #include "Offline/RecoDataProducts/inc/StrawHitFlag.hh" #include "Offline/RecoDataProducts/inc/ComboHit.hh" +#include "Offline/Mu2eKinKal/inc/KKSHFlag.hh" #include "Offline/DataProducts/inc/StrawId.hh" #include "Offline/DataProducts/inc/StrawEnd.hh" #include "Offline/DataProducts/inc/TrkTypes.hh" @@ -68,7 +69,6 @@ namespace mu2e { // calculate the doca and phi relative to the straw envelope at POCA to wire auto ppoca = XYZVectorF(uptca.particlePoca().Vect()); static XYZVectorF zdir(0.0,0.0,1.0); // relative to Z - auto wmid = XYZVectorF(straw.wirePosition()); auto wdir = XYZVectorF(straw.wireDirection()); auto delta = ppoca - wmid; @@ -87,8 +87,11 @@ namespace mu2e { if(raddir.Dot(smid) < 0.0) raddir *= -1.0; // sign radially outwards _ustrawphi = atan2(cperp.Dot(raddir),cperp.Dot(zdir)); // angle around wire WRT z axis in range -pi,pi _ustrawdist = sqrt(cperp.mag2()); - _wdot = uptca.particleDirection().Dot(straw.wireDirection()); + // set unbiased residual flags + if(udresid.active())_kkshflag.merge(KKSHFlag::goodudresid); + if(utresid.active())_kkshflag.merge(KKSHFlag::goodutresid); + if(ulresid.active())_kkshflag.merge(KKSHFlag::goodulresid); } //Legacy constructor for BTrk @@ -114,22 +117,23 @@ namespace mu2e { auto const& flag() const { return _flag; } auto const& algorithm() const { return _algo; } auto strawHitState() const { return _ambig; } - auto hitTime() const { return _etime[_eend]; } + auto time() const { return _etime[_eend]; } auto energyDep() const { return _edep; } auto const& earlyEnd() const { return _eend; } + StrawEnd lateEnd() const { return _eend.otherEnd(); } auto wireDist() const { return _wdist; } auto wireRes() const { return _werr; } auto TOTDriftTime() const { return _tottdrift; } - auto particleTOCA() const { return _ptoca; } - auto sensorTOCA() const { return _stoca; } + auto particleToca() const { return _ptoca; } + auto sensorToca() const { return _stoca; } auto fitDOCA() const { return _udoca; } auto fitDOCAVar() const { return _udocavar; } auto fitDt() const { return _udt; } - auto fitTOCAVar() const { return _utocavar; } + auto fitTocaVar() const { return _utocavar; } auto refDOCA() const { return _rdoca; } auto refDOCAVar() const { return _rdocavar; } auto refDt() const { return _rdt; } - auto reTOCAVar() const { return _rtocavar; } + auto reTocaVar() const { return _rtocavar; } auto refPOCA_Upos() const { return _rupos; } auto driftRadius() const { return _rdrift; } auto radialErr() const { return _sderr; } @@ -139,6 +143,10 @@ namespace mu2e { float signalTime() const { return _stime; } float wireDOCA() const { return _rdoca; } int ambig() const { return _ambig; } + // return a true WireHitState + WireHitState wireHitState() const { + return WireHitState(static_cast(_ambig),static_cast(_algo),_kkshflag); + } // // Payload // diff --git a/RecoDataProducts/src/KalSeed.cc b/RecoDataProducts/src/KalSeed.cc index c56befcea3..16a8432e73 100644 --- a/RecoDataProducts/src/KalSeed.cc +++ b/RecoDataProducts/src/KalSeed.cc @@ -2,13 +2,18 @@ #include namespace mu2e { + const double KalSeed::_regrowtol(1e-3); // 1 ps minimum for regrown segments + KalSeed::LHPTPtr KalSeed::loopHelixFitTrajectory() const { if(loopHelixFit() && segments().size() > 0){ // initialize the piecewise trajectory with the front segment LHPTPtr ptraj(new KalSeed::LHPT(segments().front().loopHelix())); auto iseg = segments().begin(); ++iseg; while(iseg != segments().end()){ - if(!iseg->timeRange().null())ptraj->append(iseg->loopHelix()); + if(iseg->timeRange().range() > _regrowtol ){ + auto trajptr = std::make_shared(iseg->loopHelix()); + ptraj->add(trajptr); // note this call resolves the phi0 ambiguity + } ++iseg; } return ptraj; @@ -22,7 +27,10 @@ namespace mu2e { CHPTPtr ptraj(new KalSeed::CHPT(segments().front().centralHelix())); auto iseg = segments().begin(); ++iseg; while(iseg != segments().end()){ - if(!iseg->timeRange().null())ptraj->append(iseg->centralHelix()); + if(iseg->timeRange().range() > _regrowtol ){ + auto trajptr = std::make_shared(iseg->centralHelix()); + ptraj->add(trajptr); + } ++iseg; } return ptraj; @@ -36,7 +44,10 @@ namespace mu2e { KLPTPtr ptraj(new KalSeed::KLPT(segments().front().kinematicLine())); auto iseg = segments().begin(); ++iseg; while(iseg != segments().end()){ - if(!iseg->timeRange().null())ptraj->append(iseg->kinematicLine()); + if(iseg->timeRange().range() > _regrowtol ){ + auto trajptr = std::make_shared(iseg->kinematicLine()); + ptraj->add(trajptr); + } ++iseg; } return ptraj; @@ -114,6 +125,7 @@ namespace mu2e { // start with the middle t0 = timeRange().mid(); auto iseg = nearestSegment(t0); + t0 = iseg->t0Val(_status); auto oldseg = segments().end(); unsigned ntest(0); while((!iseg->timeRange().inRange(t0)) && iseg != oldseg && ntest < segments().size()){ diff --git a/TrackerConditions/inc/DriftInfo.hh b/TrackerConditions/inc/DriftInfo.hh index 0b604b9ba0..49d4ec893c 100644 --- a/TrackerConditions/inc/DriftInfo.hh +++ b/TrackerConditions/inc/DriftInfo.hh @@ -10,7 +10,7 @@ namespace mu2e { double cDrift_ =0; // single cluster drift distance double signedDriftError_ =0; // estimated error on signed drift distance (includes LR ambiguity error effects) double unsignedDriftError_ =0; // estimated error on unsigned drift distance - double driftVelocity_; // instantaneous drift velocity + double driftVelocity_ =0; // instantaneous drift velocity static double maxdvar_; // maximum distance variance, given by straw radius double driftHitVar() const { return signedDriftError_*signedDriftError_; } // variance for hits constrained to the signed drift distance double nullHitVar() const; // variance for hits constrained to the wire position (null hits) diff --git a/TrackerConditions/test/DriftCalibPDF.C b/TrackerConditions/test/DriftCalibPDF.C index 256597b1a2..6484da00a6 100644 --- a/TrackerConditions/test/DriftCalibPDF.C +++ b/TrackerConditions/test/DriftCalibPDF.C @@ -19,277 +19,277 @@ void DriftCalibPDF(std::string filename, std::string calib_type) { #ifndef __CINT__ - TH1::AddDirectory(false); - - - // Histograms - double binsize = 0.025; - double doff = binsize * 10.; - unsigned int finefactor = 100; - - // double rdbinedges[2] = {0.0, 2.5}; - double rdbinedges[2] = {0.0, 3.0}; - int nrdbins = (rdbinedges[1] - rdbinedges[0]) / binsize; - - double rcbinedges[2] = {rdbinedges[0] - doff, rdbinedges[1] + doff}; - int nrcbins = (rcbinedges[1] - rcbinedges[0]) / binsize; - - TH1D* h_mc_dist = new TH1D("h_mc_dist", "", nrdbins, rdbinedges[0], rdbinedges[1]); // Unsigned MC true DOCA between particle trajectory and wire axis - TH1D* h_rdrift = new TH1D("h_rdrift", "", nrdbins, rdbinedges[0], rdbinedges[1]); // Cluster drift distance calibrated under old method - TH1D* h_cdrift = new TH1D("h_cdrift", "", nrcbins, rcbinedges[0], rcbinedges[1]); // Not calibrated cluster drift distance - - TH1D* h_cdrift_fine = new TH1D("h_cdrift_fine", "", nrcbins*finefactor, rcbinedges[0], rcbinedges[1]); - TH1D* h_correction = new TH1D("h_correction", "", nrcbins, rcbinedges[0], rcbinedges[1]); // Bin-by-bin corrections to not-calibrated cluster drift distance - TH1D* h_cdrift_corr = new TH1D("h_cdrift_corr", "", nrcbins, rcbinedges[0], rcbinedges[1]); // Cluster drift distance corrected - - - // Set up RooUtil - RooUtil util(filename); - - - // Event and track counters - int events_all = util.GetNEvents(); - int events_withtracks = 0; - - int tracks_all = 0; - int tracks_passcut1 = 0; - int tracks_passcut2 = 0; - int tracks_passcut3 = 0; - int tracks_goodtrackhits = 0; - - int trackhits_all = 0; - int trackhits_passcut4 = 0; - int trackhits_passcut5 = 0; - int trackhits_passcut6 = 0; - int trackhits_passcut7 = 0; - - - // Loop over events + TH1::AddDirectory(false); + + + // Histograms + double binsize = 0.025; + double doff = binsize * 10.; + unsigned int finefactor = 100; + + // double rdbinedges[2] = {0.0, 2.5}; + double rdbinedges[2] = {0.0, 3.0}; + int nrdbins = (rdbinedges[1] - rdbinedges[0]) / binsize; + + double rcbinedges[2] = {rdbinedges[0] - doff, rdbinedges[1] + doff}; + int nrcbins = (rcbinedges[1] - rcbinedges[0]) / binsize; + + TH1D* h_mc_dist = new TH1D("h_mc_dist", "", nrdbins, rdbinedges[0], rdbinedges[1]); // Unsigned MC true DOCA between particle trajectory and wire axis + TH1D* h_rdrift = new TH1D("h_rdrift", "", nrdbins, rdbinedges[0], rdbinedges[1]); // Cluster drift distance calibrated under old method + TH1D* h_cdrift = new TH1D("h_cdrift", "", nrcbins, rcbinedges[0], rcbinedges[1]); // Not calibrated cluster drift distance + + TH1D* h_cdrift_fine = new TH1D("h_cdrift_fine", "", nrcbins*finefactor, rcbinedges[0], rcbinedges[1]); + TH1D* h_correction = new TH1D("h_correction", "", nrcbins, rcbinedges[0], rcbinedges[1]); // Bin-by-bin corrections to not-calibrated cluster drift distance + TH1D* h_cdrift_corr = new TH1D("h_cdrift_corr", "", nrcbins, rcbinedges[0], rcbinedges[1]); // Cluster drift distance corrected + + + // Set up RooUtil + RooUtil util(filename); + + + // Event and track counters + int events_all = util.GetNEvents(); + int events_withtracks = 0; + + int tracks_all = 0; + int tracks_passcut1 = 0; + int tracks_passcut2 = 0; + int tracks_passcut3 = 0; + int tracks_goodtrackhits = 0; + + int trackhits_all = 0; + int trackhits_passcut4 = 0; + int trackhits_passcut5 = 0; + int trackhits_passcut6 = 0; + int trackhits_passcut7 = 0; + + + // Loop over events + // ---------------- + for ( int i_event = 0; i_event < util.GetNEvents(); ++i_event ) + { + // Get event + auto& event = util.GetEvent(i_event); + + + // Get tracks + auto tracks = event.GetTracks(is_e_minus); + + if ( event.CountTracks() <= 0 ) continue; + events_withtracks++; + + + // Loop over tracks // ---------------- - for ( int i_event = 0; i_event < util.GetNEvents(); ++i_event ) - { - // Get event - auto& event = util.GetEvent(i_event); - - - // Get tracks - auto tracks = event.GetTracks(is_e_minus); - - if ( event.CountTracks() <= 0 ) continue; - events_withtracks++; - - - // Loop over tracks - // ---------------- - for ( auto& track : tracks ) { - tracks_all++; - - // Get track particles - auto particles = track.GetMCParticles(is_track_particle); - - - // FIRST CUT: Adapted from 'gfit': ["dem.status>0"] - if ( !(track.trk->status > 0) ) continue; - tracks_passcut1++; - - // SECOND CUT: Also adapted from 'gfit' ["demmcsim[0].prirel._rel==0"] - if ( !(particles.at(0).mcsim->prirel.relationship() == 0) ) continue; - tracks_passcut2++; - - - // Get track segments - auto trksegments = track.GetSegments(); - - // Loop over track segments - bool bad_trksegment = false; - for ( auto& trksegment : trksegments ) { - if ( !has_mc_step(trksegment) ) continue; - - // SECOND CUT: Adapted from 'gmom': ["demmcsim[0].mom.R()-demmcvd.mom.R()<2.0"] - if ( !(particles.at(0).mcsim->mom.R() - trksegment.trksegmc->mom.R() < 2.0) ) { - bad_trksegment = true; - break; - } - } - if ( bad_trksegment ) continue; - tracks_passcut3++; - - - // Get track hits - auto trkhits = track.GetHits(); - - - // Loop over track hits - // -------------------- - bool passcut4 = false; - bool passcut5 = false; - bool passcut6 = false; - bool passcut7 = false; - - for ( auto& trkhit : trkhits ) { - trackhits_all++; - - // THIRD CUT: If track hit has no reconstructed object associated, skip it - if ( !(trkhit.reco != nullptr) ) continue; - bool passcut3 = true; - trackhits_passcut4++; - - // FOURTH CUT: Adapted from 'ghit': ["demtsh.state>-2"] - if ( !(trkhit.reco->state > -2) ) continue; - bool passcut4 = true; - trackhits_passcut5++; - - // FIFTH CUT: Adapted from 'thit': ["demtshmc.rel._rel==0"] - if ( !(trkhit.mc->rel.relationship() == 0) ) continue; - bool passcut5 = true; - trackhits_passcut6++; - - // SIXTH CUT: Drift quality - if ( trkhit.reco->driftqual <= 0.2 ) continue; - bool passcut6 = true; - trackhits_passcut7++; - - - // Define variables - double mc_dist = trkhit.mc->dist; - double rdrift = trkhit.reco->rdrift; - double cdrift = trkhit.reco->cdrift; - - - // Fill histograms - h_mc_dist -> Fill(mc_dist); - h_rdrift -> Fill(rdrift); - h_cdrift -> Fill(cdrift); - h_cdrift_fine -> Fill(cdrift); - - } // End of loop over track hits - - tracks_goodtrackhits++; - - } // End of loop over tracks - } // End of loop over events - - std::cout << std::endl; - std::cout << " All events: " << events_all << std::endl; - std::cout << " Events with tracks: " << events_withtracks << std::endl; - std::cout << std::endl; - std::cout << " All tracks: " << tracks_all << std::endl; - std::cout << " Tracks passing cut #1: " << tracks_passcut1 << std::endl; - std::cout << " Tracks passing cut #2: " << tracks_passcut2 << std::endl; - std::cout << " Tracks passing cut #3: " << tracks_passcut3 << std::endl; - std::cout << " Tracks with at least one good trackhit: " << tracks_goodtrackhits << std::endl; - std::cout << std::endl; - std::cout << " All track hits: " << trackhits_all << std::endl; - std::cout << " Track hits passing cut #4: " << trackhits_passcut4 << std::endl; - std::cout << " Track hits passing cut #5: " << trackhits_passcut5 << std::endl; - std::cout << " Track hits passing cut #6: " << trackhits_passcut6 << std::endl; - std::cout << " Track hits passing cut #7: " << trackhits_passcut7 << std::endl; - std::cout << std::endl; - - - // Normalize histograms - h_mc_dist -> Scale(1.0 / h_mc_dist->Integral()); - h_rdrift -> Scale(1.0 / h_rdrift->Integral()); - h_cdrift -> Scale(1.0 / h_cdrift->Integral()); - h_cdrift_fine -> Scale(1.0 / h_cdrift_fine->Integral()); - - - // Find corrections to cluster drift distance 'cdrift' - // --------------------------------------------------- - - std::vector offsets(nrcbins, 0.0); // Offset vector with length equal to number of bins of 'h_cdrift' all initialized to 0 - - int jbin = 1; - double ncalib = h_mc_dist->GetBinContent(jbin); // Content of a particular bin of 'h_mc_dist' - double cedge = h_mc_dist->GetBinLowEdge(jbin); // Lower edge of a particular bin of 'h_mc_dist' - - double calibsum = ncalib; // Sum of bin contents of 'h_mc_dist' - double nraw = 0.0; - double rawsum = nraw; // Sum of bin contents of 'h_cdrift' - - // Loop over bins of 'h_cdrift' - for ( int ibin = 1; ibin <= nrcbins; ++ibin ) { - double bincontent = h_cdrift->GetBinContent(ibin); // Bin content - double bincenter = h_cdrift->GetBinCenter(ibin); // Bin center - - rawsum += bincontent; // Add bin content to 'rawsum' - offsets[ibin-1] = bincenter - (cedge + binsize*((nraw + 0.5*bincontent)/ncalib)); // Fill offset vector entry corresponding to this bin - - // If 'nraw' plus 'h_cdrift' bin content is less than 'h_mc_dist' bin content... - if ( nraw + bincontent < ncalib ) { - nraw += bincontent; // Add 'h_cdrift' bin content to 'nraw' - } - else { - nraw += bincontent - ncalib; // Otherwise, add to 'nraw' the difference between bin contents of 'h_cdrift' minus 'h_mc_dist' - ++jbin; // Go to the next bin of 'h_mc_dist' - - // If the bin number of 'h_mc_dist' is too big, or the corresponding bin content is less than zero, break the loop - if ( jbin > nrdbins || h_mc_dist->GetBinContent(jbin) <= 0 ) break; - - // Otherwise, prepare 'ncalib' and 'cedge' with the info from the next bin of 'h_mc_dist' - ncalib = h_mc_dist->GetBinContent(jbin); - cedge = h_mc_dist->GetBinLowEdge(jbin); - - calibsum += ncalib; // Add to 'calibsum' the content of such bin of 'h_mc_dist' + for ( auto& track : tracks ) { + tracks_all++; + + // Get track particles + auto particles = track.GetMCParticles(is_track_particle); + + + // FIRST CUT: Adapted from 'gfit': ["dem.status>0"] + if ( !(track.trk->status > 0) ) continue; + tracks_passcut1++; + + // SECOND CUT: Also adapted from 'gfit' ["demmcsim[0].prirel._rel==0"] + if ( !(particles.at(0).mcsim->prirel.relationship() == 0) ) continue; + tracks_passcut2++; + + + // Get track segments + auto trksegments = track.GetSegments(); + + // Loop over track segments + bool bad_trksegment = false; + for ( auto& trksegment : trksegments ) { + if ( !has_mc_step(trksegment) ) continue; + + // SECOND CUT: Adapted from 'gmom': ["demmcsim[0].mom.R()-demmcvd.mom.R()<2.0"] + if ( !(particles.at(0).mcsim->mom.R() - trksegment.trksegmc->mom.R() < 2.0) ) { + bad_trksegment = true; + break; } + } + if ( bad_trksegment ) continue; + tracks_passcut3++; + + + // Get track hits + auto trkhits = track.GetHits(); + + + // Loop over track hits + // -------------------- + bool passcut4 = false; + bool passcut5 = false; + bool passcut6 = false; + bool passcut7 = false; + + for ( auto& trkhit : trkhits ) { + trackhits_all++; + + // THIRD CUT: If track hit has no reconstructed object associated, skip it + if ( !(trkhit.reco != nullptr) ) continue; + bool passcut3 = true; + trackhits_passcut4++; + + // FOURTH CUT: Adapted from 'ghit': ["demtsh.state>-2"] + if ( !(trkhit.reco->state > -2) ) continue; + bool passcut4 = true; + trackhits_passcut5++; + + // FIFTH CUT: Adapted from 'thit': ["demtshmc.rel._rel==0"] + if ( !(trkhit.mc->rel.relationship() == 0) ) continue; + bool passcut5 = true; + trackhits_passcut6++; + + // SIXTH CUT: Drift quality + if ( trkhit.reco->driftqual <= 0.2 ) continue; + bool passcut6 = true; + trackhits_passcut7++; + + + // Define variables + double mc_dist = trkhit.mc->dist; + double rdrift = trkhit.reco->rdrift; + double cdrift = trkhit.reco->cdrift; + + + // Fill histograms + h_mc_dist -> Fill(mc_dist); + h_rdrift -> Fill(rdrift); + h_cdrift -> Fill(cdrift); + h_cdrift_fine -> Fill(cdrift); + + } // End of loop over track hits + + tracks_goodtrackhits++; + + } // End of loop over tracks + } // End of loop over events + + std::cout << std::endl; + std::cout << " All events: " << events_all << std::endl; + std::cout << " Events with tracks: " << events_withtracks << std::endl; + std::cout << std::endl; + std::cout << " All tracks: " << tracks_all << std::endl; + std::cout << " Tracks passing cut #1: " << tracks_passcut1 << std::endl; + std::cout << " Tracks passing cut #2: " << tracks_passcut2 << std::endl; + std::cout << " Tracks passing cut #3: " << tracks_passcut3 << std::endl; + std::cout << " Tracks with at least one good trackhit: " << tracks_goodtrackhits << std::endl; + std::cout << std::endl; + std::cout << " All track hits: " << trackhits_all << std::endl; + std::cout << " Track hits passing cut #4: " << trackhits_passcut4 << std::endl; + std::cout << " Track hits passing cut #5: " << trackhits_passcut5 << std::endl; + std::cout << " Track hits passing cut #6: " << trackhits_passcut6 << std::endl; + std::cout << " Track hits passing cut #7: " << trackhits_passcut7 << std::endl; + std::cout << std::endl; + + + // Normalize histograms + h_mc_dist -> Scale(1.0 / h_mc_dist->Integral()); + h_rdrift -> Scale(1.0 / h_rdrift->Integral()); + h_cdrift -> Scale(1.0 / h_cdrift->Integral()); + h_cdrift_fine -> Scale(1.0 / h_cdrift_fine->Integral()); + + + // Find corrections to cluster drift distance 'cdrift' + // --------------------------------------------------- + + std::vector offsets(nrcbins, 0.0); // Offset vector with length equal to number of bins of 'h_cdrift' all initialized to 0 + + int jbin = 1; + double ncalib = h_mc_dist->GetBinContent(jbin); // Content of a particular bin of 'h_mc_dist' + double cedge = h_mc_dist->GetBinLowEdge(jbin); // Lower edge of a particular bin of 'h_mc_dist' + + double calibsum = ncalib; // Sum of bin contents of 'h_mc_dist' + double nraw = 0.0; + double rawsum = nraw; // Sum of bin contents of 'h_cdrift' + + // Loop over bins of 'h_cdrift' + for ( int ibin = 1; ibin <= nrcbins; ++ibin ) { + double bincontent = h_cdrift->GetBinContent(ibin); // Bin content + double bincenter = h_cdrift->GetBinCenter(ibin); // Bin center + + rawsum += bincontent; // Add bin content to 'rawsum' + offsets[ibin-1] = bincenter - (cedge + binsize*((nraw + 0.5*bincontent)/ncalib)); // Fill offset vector entry corresponding to this bin + + // If 'nraw' plus 'h_cdrift' bin content is less than 'h_mc_dist' bin content... + if ( nraw + bincontent < ncalib ) { + nraw += bincontent; // Add 'h_cdrift' bin content to 'nraw' } - - - // Perform a safety check - if ( jbin < nrcbins && (std::fabs(1.0 - calibsum) > 1.0e-6 || std::fabs(1.0 - rawsum) > 1.0e-6) ) { - std::cout << " Early exit! " << std::endl; - std::cout << " \tjbin = " << jbin << std::endl; - std::cout << " \tRemaining MC fraction = " << 1.0 - calibsum << std::endl; - std::cout << " \tcdrift fraction = " << 1.0 - calibsum << std::endl; - } - - - // Fill bin-by-bin corrections and corrected drift histograms - for ( unsigned int ibin = 0; ibin < offsets.size(); ++ibin ) { - h_correction -> SetBinContent(ibin+1, offsets[ibin]); - } - - for ( unsigned int ibin = 1; ibin <= finefactor*nrcbins; ++ibin ) { - unsigned int jbin = floor(float(ibin)/float(finefactor)); - h_cdrift_corr -> Fill(h_cdrift_fine->GetBinCenter(ibin+1) - offsets[jbin], h_cdrift_fine->GetBinContent(ibin+1)); - } - h_cdrift_corr -> Scale(1.0 / h_cdrift_corr->Integral()); - - - // Save calibration file - std::string cfname = "DriftCalibPDF"; - ofstream cfile(Form("%s_%s.txt", cfname.c_str(), calib_type.c_str()), ios::trunc); - time_t now = time(0); - char* dt = ctime(&now); - - cfile << "# The following was produced by DriftCalibPDF.C with track hit selection 'trkhit.reco->driftqual > 0.2' on " << dt << std::endl; - cfile << std::setw(4) << std::setprecision(3); - cfile << "driftOffBins : [ "; - cfile << rcbinedges[0] << " , " << rcbinedges[1] << " ]" << endl; - cfile << "driftOffset : [ "; - - bool first(true); - for ( int ibin = 0; ibin < nrcbins; ++ibin ) { - if ( !first ) cfile << " , "; - first = false; - cfile << offsets[ibin]; + else { + nraw += bincontent - ncalib; // Otherwise, add to 'nraw' the difference between bin contents of 'h_cdrift' minus 'h_mc_dist' + ++jbin; // Go to the next bin of 'h_mc_dist' + + // If the bin number of 'h_mc_dist' is too big, or the corresponding bin content is less than zero, break the loop + if ( jbin > nrdbins || h_mc_dist->GetBinContent(jbin) <= 0 ) break; + + // Otherwise, prepare 'ncalib' and 'cedge' with the info from the next bin of 'h_mc_dist' + ncalib = h_mc_dist->GetBinContent(jbin); + cedge = h_mc_dist->GetBinLowEdge(jbin); + + calibsum += ncalib; // Add to 'calibsum' the content of such bin of 'h_mc_dist' } - cfile << " ]" << endl; - cfile.close(); - - - // Save histograms - TFile outputfile(Form("../root_files/DriftCalibPDF_%s.root", calib_type.c_str()), "RECREATE"); - - h_mc_dist -> Write(); - h_rdrift -> Write(); - h_cdrift -> Write(); - - h_cdrift_fine -> Write(); - h_correction -> Write(); - h_cdrift_corr -> Write(); - - outputfile.Close(); + } + + + // Perform a safety check + if ( jbin < nrcbins && (std::fabs(1.0 - calibsum) > 1.0e-6 || std::fabs(1.0 - rawsum) > 1.0e-6) ) { + std::cout << " Early exit! " << std::endl; + std::cout << " \tjbin = " << jbin << std::endl; + std::cout << " \tRemaining MC fraction = " << 1.0 - calibsum << std::endl; + std::cout << " \tcdrift fraction = " << 1.0 - calibsum << std::endl; + } + + + // Fill bin-by-bin corrections and corrected drift histograms + for ( unsigned int ibin = 0; ibin < offsets.size(); ++ibin ) { + h_correction -> SetBinContent(ibin+1, offsets[ibin]); + } + + for ( unsigned int ibin = 1; ibin <= finefactor*nrcbins; ++ibin ) { + unsigned int jbin = floor(float(ibin)/float(finefactor)); + h_cdrift_corr -> Fill(h_cdrift_fine->GetBinCenter(ibin+1) - offsets[jbin], h_cdrift_fine->GetBinContent(ibin+1)); + } + h_cdrift_corr -> Scale(1.0 / h_cdrift_corr->Integral()); + + + // Save calibration file + std::string cfname = "DriftCalibPDF"; + ofstream cfile(Form("%s_%s.txt", cfname.c_str(), calib_type.c_str()), ios::trunc); + time_t now = time(0); + char* dt = ctime(&now); + + cfile << "# The following was produced by DriftCalibPDF.C with track hit selection 'trkhit.reco->driftqual > 0.2' on " << dt << std::endl; + cfile << std::setw(4) << std::setprecision(3); + cfile << "driftOffBins : [ "; + cfile << rcbinedges[0] << " , " << rcbinedges[1] << " ]" << endl; + cfile << "driftOffset : [ "; + + bool first(true); + for ( int ibin = 0; ibin < nrcbins; ++ibin ) { + if ( !first ) cfile << " , "; + first = false; + cfile << offsets[ibin]; + } + cfile << " ]" << endl; + cfile.close(); + + + // Save histograms + TFile outputfile(Form("../root_files/DriftCalibPDF_%s.root", calib_type.c_str()), "RECREATE"); + + h_mc_dist -> Write(); + h_rdrift -> Write(); + h_cdrift -> Write(); + + h_cdrift_fine -> Write(); + h_correction -> Write(); + h_cdrift_corr -> Write(); + + outputfile.Close(); #endif // __CINT__ } diff --git a/TrkFilters/src/KalSeedFilter_module.cc b/TrkFilters/src/KalSeedFilter_module.cc index 6dd2d6ffd4..acaa3edb95 100644 --- a/TrkFilters/src/KalSeedFilter_module.cc +++ b/TrkFilters/src/KalSeedFilter_module.cc @@ -182,10 +182,10 @@ namespace mu2e bool KalSeedFilter::checkKalSeed(const KalSeed&Ks, const KalSeedCutsTool&Cuts){ if(_debug > 3){ - std::cout << "KalSeedFilter: in checkKalSeed status "<< Ks.status() << " intersections " << Ks.intersections().size() << std::endl; + std::cout << "KalSeedFilter: in checkKalSeed status "<< Ks.status() << std::endl; } - if( Ks.status().hasAllProperties(Cuts._goods) && Ks.intersections().size()>0){ + if( Ks.status().hasAllProperties(Cuts._goods) ){ // extract test quantities from the fit segment at t0 double t0; diff --git a/TrkReco/src/TrackMatching_module.cc b/TrkReco/src/TrackMatching_module.cc index b8ab892c7c..8cd47ee53e 100644 --- a/TrkReco/src/TrackMatching_module.cc +++ b/TrkReco/src/TrackMatching_module.cc @@ -109,7 +109,7 @@ namespace mu2e if(_useIndices) { // match using the digi index if(hit_1._index == hit_2._index) ++overlap; } else { - if(hit_1._sid == hit_2._sid && std::fabs(hit_1.hitTime() - hit_2.hitTime()) < _timeWindow) ++overlap; + if(hit_1._sid == hit_2._sid && std::fabs(hit_1.time() - hit_2.time()) < _timeWindow) ++overlap; } } }