diff --git a/DataProducts/inc/SurfaceId.hh b/DataProducts/inc/SurfaceId.hh index 6147a3150e..860faed1b5 100644 --- a/DataProducts/inc/SurfaceId.hh +++ b/DataProducts/inc/SurfaceId.hh @@ -38,12 +38,10 @@ namespace mu2e { CRV_StrongBack=280, // CRV module Al strongback (tracker-side support plate) DS_HatchConcrete=300 // detector-area hatch concrete block approximation }; - - // Update this counter whenever you add/remove surface IDs from the enum above. - static constexpr std::size_t nSurfaceIds = 64; - - static std::string const& typeName(); - static std::map const& names(); + // Update this counter whenever you add/remove surface IDs from the enum above. + static constexpr std::size_t nSurfaceIds = 64; + static std::string const& typeName(); + static std::map const& names(); }; using SurfaceIdEnum = EnumToStringSparse; class SurfaceId { @@ -59,12 +57,13 @@ namespace mu2e { int index() const { return index_; } auto const& name() const { return sid_.name(); } - bool indexMatch(SurfaceId const& other) const { return index_ == other.index_ || index_ < 0 || other.index_ < 0; } - bool indexCompare(SurfaceId const& other) const { return index_<0 || other.index_ < 0 ? false : index_ < other.index_; } + bool indexMatch(SurfaceId const& other) const { return index_ == other.index_ || index_ == allIndices_ || other.index_ == allIndices_; } + bool indexCompare(SurfaceId const& other) const { return index_== allIndices_ || other.index_ == allIndices_ ? false : index_ < other.index_; } bool operator == (SurfaceId const& other ) const { return sid_ == other.sid_ && indexMatch(other) ; } bool operator != (SurfaceId const& other ) const { return sid_ != other.sid_ || !indexMatch(other) ; } bool operator < (SurfaceId const& other ) const { return sid_ == other.sid_ ? indexCompare(other) : sid_ < other.sid_; } - private: + static int allIndices_; // index wildcard + private: SurfaceIdEnum sid_; int index_; // index. Negative value is a wild card for matching }; diff --git a/DataProducts/src/SurfaceId.cc b/DataProducts/src/SurfaceId.cc index bef10f3b4f..e4bdfe2355 100644 --- a/DataProducts/src/SurfaceId.cc +++ b/DataProducts/src/SurfaceId.cc @@ -12,6 +12,7 @@ namespace mu2e { static const std::string type("SurfaceIdEnum"); return type; } + int SurfaceId::allIndices_(-1); // wildcard for index comparison namespace { using SurfaceIdName = std::pair; diff --git a/KinKalGeom/CMakeLists.txt b/KinKalGeom/CMakeLists.txt index b99bed90b4..f9014a79d7 100644 --- a/KinKalGeom/CMakeLists.txt +++ b/KinKalGeom/CMakeLists.txt @@ -5,6 +5,7 @@ cet_make_library( src/KKStrawMaterial.cc LIBRARIES PUBLIC Offline::ConfigTools + Offline::DataProducts KinKal::MatEnv KinKal::Geometry KinKal::Trajectory diff --git a/KinKalGeom/fcl/prolog.fcl b/KinKalGeom/fcl/prolog.fcl index 8af7a6f49e..2a15bbdfcc 100644 --- a/KinKalGeom/fcl/prolog.fcl +++ b/KinKalGeom/fcl/prolog.fcl @@ -13,7 +13,7 @@ KinKalGeom : { IonizationEnergyLossMode : 1 # Moyal mean SolidScatteringFraction : 0.999999 GasScatteringFraction : 0.9999999 - ElectronBrehmsFraction : 0.04 + ElectronBrehmsFraction : 0.03 # Bethe path correction for KKShellXing: off by default; enable per crossing type as needed. # Only effective when IonizationEnergyLossMode is moyalmean. BetheCorrectionIPA : false diff --git a/KinKalGeom/src/SConscript b/KinKalGeom/src/SConscript index eea1886bbb..8144b196ad 100644 --- a/KinKalGeom/src/SConscript +++ b/KinKalGeom/src/SConscript @@ -15,6 +15,7 @@ helper = mu2e_helper(env) mainlib = helper.make_mainlib([ 'mu2e_ConfigTools', + 'mu2e_DataProducts', 'KinKal_MatEnv', 'KinKal_Geometry', 'KinKal_Trajectory', diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 9d4acebcc1..42c9b45cfc 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -636,9 +636,10 @@ Mu2eKinKal : { } ExtrapolationSettings :{ @table::Mu2eKinKal.LHDRIFTXTRAP - BackToTracker: true + BackToTracker : true } Extend: false + MustRegrow : "Any" } producers : { diff --git a/Mu2eKinKal/inc/ExtrapolateIPA.hh b/Mu2eKinKal/inc/ExtrapolateIPA.hh index 0b583ba505..59eaec42f0 100644 --- a/Mu2eKinKal/inc/ExtrapolateIPA.hh +++ b/Mu2eKinKal/inc/ExtrapolateIPA.hh @@ -31,6 +31,7 @@ namespace mu2e { auto const& intersection() const { return inter_; } double zmin() const { return zmin_; } double zmax() const { return zmax_; } + double zmid() const { return 0.5*(zmin_+zmax_); } int debug() const { return debug_; } // extrapolation predicate: the track will be extrapolated until this predicate returns false, subject to the maximum time template bool needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const; @@ -50,47 +51,41 @@ namespace mu2e { }; template bool ExtrapolateIPA::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const { - // we are answering the question: did the segment last added to this extrapolated track hit the IPA or not? - // if so, stop extrapolating (for now). If not, and if we're still inside or heading towards the IPA, keep going. + // we are answering the question: did the segment last added to this extrapolated trajectory hit the IPA or not? + // if so, stop extrapolating (for now). If not, and if we're still inside or heading towards the IPA, keep going, otherwise stop reset(); // clear any cache - auto const& ktraj = tdir == TimeDir::forwards ? fittraj.back() : fittraj.front(); + bool retval(true); // by default keep going + auto const& ktraj = tdir == TimeDir::forwards ? fittraj.back() : fittraj.front(); // most recently added segment // add a small buffer to the test range to prevent re-intersection with the same piece static const double epsilon(1e-7); // small step to avoid re-intersecting 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.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 - if(debug_ > 2)std::cout << "IPA extrap start time " << stime << " start z " << spos.Z() << " end z " << epos.Z() << " zvel " << zvel << std::endl; - // stop if the particle is heading away from the IPA - if( (zvel > 0 && spos.Z() > zmax_ ) || (zvel < 0 && spos.Z() < zmin_)){ - if(debug_ > 1)std::cout << "Heading away from IPA: done" << std::endl; - return false; - } - // if the particle is going in the right direction but hasn't yet reached the IPA just keep going - if( (zvel > 0 && epos.Z() < zmin_) || (zvel < 0 && epos.Z() > zmax_) ){ - if(debug_ > 2)std::cout << "Heading towards IPA, z " << spos.Z()<< std::endl; - return true; - } - // if we get to here we need to test for an intersection with the actual cylinder. Make sure the range is positive definite - auto trange = tdir == TimeDir::forwards ? TimeRange(stime,etime) : TimeRange(etime,stime); - Intersection newinter = KinKal::intersect(fittraj,*ipa_,trange,intertol_,tdir); - if(debug_ > 2)std::cout << "IPA " << newinter << std::endl; - if(newinter.good()){ - // update the cache - inter_ = newinter; - if(debug_ > 0)std::cout << "Good IPA " << newinter << std::endl; - return false; + auto trange = tdir == TimeDir::forwards ? + TimeRange(ktraj.range().begin()+epsilon, ktraj.range().end()) : + TimeRange(ktraj.range().begin(), ktraj.range().end()-epsilon); + auto spos = ktraj.position3(trange.begin()); + auto epos = ktraj.position3(trange.end()); + // if either end of the segment is inside the z range of the IPA, test for an in-range interestection + if( (spos.Z() > zmin_ && spos.Z() < zmax_) || (epos.Z() > zmin_ && epos.Z() < zmax_) ){ + Intersection newinter = KinKal::intersect(ktraj,*ipa_,trange,intertol_,tdir); + if(debug_ > 2)std::cout << "IPA " << newinter << std::endl; + if(newinter.good()){ + // update the cache + inter_ = newinter; + if(debug_ > 0)std::cout << "Good IPA " << newinter << ", Stopping " << std::endl; + retval = false; + } } else { - // no more intersections: keep extending in Z till we clear the IPA - if(debug_ > 1)std::cout << "Extrapolating to IPA edge, z " << spos.Z() << std::endl; - if(zvel > 0.0) - return spos.Z() < zmax_; - else - return spos.Z() > zmin_; + // otherwise, if the trajectory is heading towards the IPA keep going + auto vel = ktraj.velocity(trange.end())*timeDirSign(tdir); // sign velocity by extrapolation direction + double vb = vel.Dot(ktraj.bnom().Unit()); // project along the BField axis. This is insensitive to reflection. + if( (epos.Z() > zmax_ && vb < 0) || (epos.Z() < zmin_ && vb > 0) ) { + if(debug_ > 1)std::cout << "Extrapolating towards IPA, z " << epos.Z() << " bvel " << vb << std::endl; + } else { + retval = false; + if(debug_ > 1)std::cout << "Heading away from IPA, z " << epos.Z() << " bvel " << vb << std::endl; + } } + return retval; } } #endif diff --git a/Mu2eKinKal/inc/ExtrapolateST.hh b/Mu2eKinKal/inc/ExtrapolateST.hh index 54341af354..8ed09a8808 100644 --- a/Mu2eKinKal/inc/ExtrapolateST.hh +++ b/Mu2eKinKal/inc/ExtrapolateST.hh @@ -36,6 +36,7 @@ namespace mu2e { double interTolerance() const { return intertol_; } double zmin() const { return zmin_; } double zmax() const { return zmax_; } + double zmid() const { return 0.5*(zmin_+zmax_); } double rmin() const { return rmin_; } double rmax() const { return rmax_; } auto const& foils () const { return foils_; } @@ -47,8 +48,8 @@ namespace mu2e { template bool needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const; // reset between tracks void reset() const { inter_ = Intersection(); sid_ = SurfaceId(); ann_ = AnnPtr();} - // find the nearest foil to a z positionin a given z direction - size_t nearestFoil(double zpos, double zdir) const; + // find the foils bounded by a set of z positions in a given z direction + void boundedFoils(double spos, double epos, double zdir, std::vector& foils) const; private: double maxDt_ = -1; // maximum extrapolation time double maxDtStep_ = -1; // maximum extrapolation time step in a single iteration @@ -70,81 +71,77 @@ namespace mu2e { // we are answering the question: did the segment last added to this extrapolated track hit a foil or not? If so we are done // extrapolating (for now) and we want to find all the intersections in that piece. If not, and if we're still inside or heading towards the // ST, keep going. + reset(); // clear any cache + bool retval = true; // by default keep going auto const& ktraj = tdir == TimeDir::forwards ? fittraj.back() : fittraj.front(); // add a small buffer to the test range to prevent re-intersection with the same piece static const double epsilon(1e-7); // small difference to avoid re-intersecting 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.velocity(stime); // physical velocity - if(tdir == TimeDir::backwards) vel *= -1.0; - auto spos = ktraj.position3(stime); - auto epos = ktraj.position3(etime); + auto trange = tdir == TimeDir::forwards ? + TimeRange(ktraj.range().begin()+epsilon, ktraj.range().end()) : + TimeRange(ktraj.range().begin(), ktraj.range().end()-epsilon); + auto spos = ktraj.position3(trange.begin()); + auto epos = ktraj.position3(trange.end()); + auto vel = ktraj.velocity(trange.begin())*timeDirSign(tdir); // sign velocity by extrapolation direction if(debug_ > 2)std::cout << "ST extrap tdir " << tdir << " start z " << spos.Z() << " end z " << epos.Z() << " zvel " << vel.Z() << " rho " << spos.Rho() << std::endl; - // stop if the particle is heading away from the ST - if( (vel.Z() > 0 && spos.Z() > zmax_ ) || (vel.Z() < 0 && spos.Z() < zmin_)){ - reset(); // clear any cache - if(debug_ > 1)std::cout << "Heading away from ST: done" << std::endl; - return false; - } - // if the particle is going in the right direction but haven't yet reached the ST in Z just keep going - if( (vel.Z() > 0 && epos.Z() < zmin_) || (vel.Z() < 0 && epos.Z() > zmax_) ){ - reset(); - if(debug_ > 2)std::cout << "Heading towards ST, z " << spos.Z()<< std::endl; - return true; - } - // if we get to here we are in the correct Z range. Test foils. - int ifoil = nearestFoil(spos.Z(),vel.Z()); - if(ifoil >= (int)foils_.size())return true; - if(debug_ > 2)std::cout << "Looping on foils " << std::endl; - int dfoil = vel.Z() > 0.0 ? 1 : -1; // iteration direction - auto trange = tdir == TimeDir::forwards ? TimeRange(stime,ktraj.range().end()) : TimeRange(ktraj.range().begin(),stime); - // loop over foils in the z range of this piece - while(ifoil >= 0 && ifoil < (int)foils_.size() && (foils_[ifoil]->center().Z() - epos.Z())*dfoil < 0.0){ - auto foilptr = foils_[ifoil]; - if(debug_ > 2)std::cout << "foil " << ifoil << " z " << foilptr->center().Z() << std::endl; - auto newinter = KinKal::intersect(ktraj,*foilptr,trange,intertol_,tdir); - if(debug_ > 2)std::cout << "ST foil inter " << newinter << std::endl; - if(newinter.good()){ - // update the cache - inter_ = newinter; - ann_ = foils_[ifoil]; - sid_ = SurfaceId(SurfaceIdEnum::ST_Foils,ifoil); - if(debug_ > 0)std::cout << "Good ST foil " << newinter << " sid " << sid_ << std::endl; - return false; + // rough geometric test + if( (spos.Z() > zmin_ && spos.Z() < zmax_) || (epos.Z() > zmin_ && epos.Z() < zmax_) ) { + double srho = spos.Rho(); + double erho = epos.Rho(); + if( ( srho > rmin_ && srho < rmax_) || (erho > rmin_ && erho < rmax_)){ + // Find foils bounded by the segment Z positions. These are ordered by the trajectory direction, so they can be tested in order + std::vector foils; + boundedFoils(spos.Z(),epos.Z(),vel.Z(),foils); + if(debug_ > 2)std::cout << "Looping on " << foils.size() << " foils " << std::endl; + // loop over foils in the z range of this piece, and find the first intersection by this piece. + for(auto ifoil : foils){ + auto const& foilptr = foils_[ifoil]; + if(debug_ > 2)std::cout << "foil " << ifoil << " z " << foilptr->center().Z() << std::endl; + auto newinter = KinKal::intersect(ktraj,*foilptr,trange,intertol_,tdir); + if(debug_ > 2)std::cout << "ST foil inter " << newinter << std::endl; + if(newinter.good()){ + // update the cache + inter_ = newinter; + ann_ = foils_[ifoil]; + sid_ = SurfaceId(SurfaceIdEnum::ST_Foils,ifoil); + if(debug_ > 0)std::cout << "Good ST foil " << newinter << " sid " << sid_ << std::endl; + retval = false; + break; + } + } + } + } else { + // if the particle is going in the right direction but hasn't yet reached the ST in Z just keep going + double vb = vel.Dot(ktraj.bnom().Unit()); + if( (epos.Z() > zmax_ && vb < 0) || (epos.Z() < zmin_ && vb > 0) ) { + if(debug_ > 1)std::cout << "Extrapolating towards ST, z " << epos.Z() << " bvel " << vb << std::endl; + } else { + retval = false; + if(debug_ > 1)std::cout << "Heading away from ST, z " << epos.Z() << " bvel " << vb << std::endl; } - ifoil += dfoil; // otherwise continue loopin on foils } - // no more intersections: keep extending in Z till we clear the ST - reset(); - if(debug_ > 1)std::cout << "Extrapolating to ST edge, z " << spos.Z() << std::endl; - if(vel.Z() > 0.0) - return spos.Z() < zmax_; - else - return spos.Z() > zmin_; + return retval; } - size_t ExtrapolateST::nearestFoil(double zpos, double zvel) const { - size_t retval = foils_.size(); + void ExtrapolateST::boundedFoils(double sz, double ez, double zvel, std::vector& foils) const { + double zmin = min(sz,ez); + double zmax = max(sz,ez); if(zvel > 0.0){ // going forwards in z - for(auto ifoil= foils_.begin(); ifoil != foils_.end(); ifoil++){ - auto const& foilptr = *ifoil; - if(foilptr->center().Z() > zpos){ - retval = std::distance(foils_.begin(),ifoil); - break; + for(size_t ifoil=0;ifoil < foils_.size();++ifoil) { + auto const& foilptr = foils_[ifoil]; + if(foilptr->center().Z() > zmin && foilptr->center().Z() < zmax ){ + foils.push_back(ifoil); } } - } else { - for(auto ifoil= foils_.rbegin(); ifoil != foils_.rend(); ifoil++){ - auto const& foilptr = *ifoil; - if(foilptr->center().Z() < zpos){ - auto jfoil = ifoil.base()-1; // points to the equivalent forwards object - retval = std::distance(foils_.begin(),jfoil); - break; + } else { // backwards in z + for(int ifoil = foils_.size()-1;ifoil >= 0; --ifoil){ + auto jfoil = static_cast(ifoil); + auto const& foilptr = foils_[jfoil]; + if(foilptr->center().Z() > zmin && foilptr->center().Z() < zmax ){ + foils.push_back(jfoil); } } } - return retval; } } diff --git a/Mu2eKinKal/inc/ExtrapolateToZ.hh b/Mu2eKinKal/inc/ExtrapolateToZ.hh index f4a4fc513f..7150e8fc1c 100644 --- a/Mu2eKinKal/inc/ExtrapolateToZ.hh +++ b/Mu2eKinKal/inc/ExtrapolateToZ.hh @@ -4,6 +4,7 @@ #define Mu2eKinKal_ExtrapolateToZ_hh #include "KinKal/Trajectory/ParticleTrajectory.hh" #include "KinKal/General/TimeDir.hh" +#include "KinKal/Geometry/Plane.hh" #include #include "cetlib_except/exception.h" namespace mu2e { @@ -11,7 +12,8 @@ namespace mu2e { using KinKal::timeDirSign; class ExtrapolateToZ { public: - ExtrapolateToZ(double maxdt, double maxdtstep, double tol, double zval,int debug) : maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(tol), zval_(zval), debug_(debug) {} + ExtrapolateToZ(double maxdt, double maxdtstep, double tol, double intertol, double zval,int debug) : maxDt_(maxdt), maxDtStep_(maxdtstep), + dptol_(tol), intertol_(intertol), zval_(zval), debug_(debug), plane_(VEC3(0.0,0.0,1.0), VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,zval_)) {} // interface for extrapolation double maxDt() const { return maxDt_; } // maximum time to extend the track, WRT the time of the first(last) measurement double maxDtStep() const { return maxDtStep_; } @@ -24,30 +26,31 @@ namespace mu2e { double maxDt_ = -1; // maximum extrapolation time double maxDtStep_ = -1; // maximum extrapolation time step in a single iteration double dptol_ = 1e10; // fractional momentum tolerance in BField domain + double intertol_ = 1e-6; // intersection distance tolerance double zval_ = 0; // z value targeted int debug_ = 0; // debug level + KinKal::Plane plane_; // KinKal plane }; template bool ExtrapolateToZ::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, KinKal::TimeDir tdir) const { auto const& ktraj = tdir == TimeDir::forwards ? fittraj.back() : fittraj.front(); auto time = tdir == TimeDir::forwards ? ktraj.range().end() : ktraj.range().begin(); - auto vel = ktraj.velocity(time); - auto pos = ktraj.position3(time); - double zvel = vel.Z()*timeDirSign(tdir); // sign by extrapolation direction - double zval = pos.Z(); - auto const& bnom = ktraj.bnom(time); - double zref = vel.R()*fabs(sin(bnom.Theta())); - if(debug_ > 2)std::cout << "Z extrap start time " << time << " z " << zval << " zvel " << zvel << " zref " << zref << std::endl; - // if z velocity is unreliable, continue - if(fabs(zvel) < zref) return true; - // stop if we're heading away from the target z - if((zvel > 0 && zval > zval_ ) || (zvel < 0 && zval < zval_))return false; - // stop when we get beyond the target value in Z - if(zvel< 0){ // backwards extrapolation of downstream-going track, or forwards extrapolation of upstream-going track - return zval > zval_; - } else { // opposite - return zval < zval_; - } + auto vel = ktraj.velocity(time)*timeDirSign(tdir); // sign by extrapolation direction + auto epos =ktraj.position3(time); + double dz = zval_ - epos.Z(); + if(debug_ > 2)std::cout << "Z extrap start time " << time << " dz " << dz << " zvel " << vel.Z() << " dz " << dz << std::endl; + // false if we're heading away from the target z, otherwise true + bool retval (true); + if(dz*vel.Z() < 0) { + // near reflection vz becomes an unreliable estimator of direction. Use the linearized direction instead + auto axis = ktraj.linearize(time); + double dist; + plane_.intersect(axis,dist,tdir==TimeDir::forwards,intertol_); + retval = dist*timeDirSign(tdir) > 0.0; + if(debug_ > 2)std::cout << "Z extrap ray intersection result " << retval << " distance " << dist << std::endl; + } else if (debug_ > 2) + std::cout << "Z heading towards surface, distance " << dz << std::endl; + return retval; } } #endif diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index 907db164e1..13bee1d9db 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -46,7 +46,7 @@ namespace mu2e { private: int debug_; double btol_, intertol_, maxdt_, maxdtstep_, minv_; - bool backToTracker_, extrapolateOPA_, toTrackerEnds_, upstream_, toCRV_; + bool backToTracker_, extrapolateOPA_, toTrackerEnds_, upstream_, toCRV_, stopAtTSDA_; 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 }; @@ -81,7 +81,7 @@ namespace mu2e { 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 + 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); @@ -95,8 +95,7 @@ namespace mu2e { } } } else { // reflection inside the IPA; extrapolate back through the IPA, then to the tracker entrance - ExtrapolateToZ trackerFront(maxdt_,maxdtstep_,btol_,kkg_h->tracker()->front().center().Z(),debug_); - if(backToTracker_)ktrk.extrapolate(tdir,trackerFront); + if(backToTracker_)extrapolateTracker(ktrk,tdir); } // optionally test for intersection with the OPA if(extrapolateOPA_)extrapolateOPA(ktrk,starttime,tdir); @@ -111,8 +110,8 @@ namespace mu2e { template void KKExtrap::toTrackerEnds(KKTrack& ktrk) const { GeomHandle kkg_h; - ExtrapolateToZ trackerFront(maxdt_,maxdtstep_,btol_,kkg_h->tracker()->front().center().Z(),debug_); - ExtrapolateToZ trackerBack(maxdt_,maxdtstep_,btol_,kkg_h->tracker()->back().center().Z(),debug_); + ExtrapolateToZ trackerFront(maxdt_,maxdtstep_,btol_,intertol_,kkg_h->tracker()->front().center().Z(),debug_); + ExtrapolateToZ trackerBack(maxdt_,maxdtstep_,btol_,intertol_,kkg_h->tracker()->back().center().Z(),debug_); // time direction to reach the bounding surfaces from the active region depends on the z momentum. This calculation assumes the particle doesn't // reflect inside the tracker volume auto const& ftraj = ktrk.fitTraj(); @@ -163,7 +162,7 @@ namespace mu2e { 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); + double startz = ftraj.position3(starttime).Z() - extrapIPA.zmid(); // z relative to the IPA center do { ktrk.extrapolate(tdir,extrapIPA); if(extrapIPA.intersection().good()){ @@ -184,13 +183,16 @@ namespace mu2e { } } } while(extrapIPA.intersection().good()); - // check if the particle exited in the same physical direction or not (reflection) + // check if the particle exited the same side it entered (ie it reflected inside the IPA) double endtime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); - auto enddir = ftraj.direction(endtime); - if(enddir.Z() * startdir.Z() > 0.0){ + double endz = ftraj.position3(endtime).Z() - extrapIPA.zmid(); // z relative to the IPA center + if(endz * startz < 0.0){ + if(extrapIPA.debug() > 1)std::cout <<"Exiting opposite side IPA" << std::endl; return true; + } else { + if(extrapIPA.debug() > 1)std::cout <<"Exiting same side IPA" << std::endl; + return false; } - return false; } template bool KKExtrap::extrapolateST(KKTrack& ktrk,TimeDir tdir) const { @@ -204,7 +206,7 @@ namespace mu2e { ExtrapolateST extrapST(maxdt_,maxdtstep_,btol_,intertol_,*kkg_h->ST(),debug_); auto const& ftraj = ktrk.fitTraj(); double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); - auto startdir = ftraj.direction(starttime); + double startz = ftraj.position3(starttime).Z() - extrapST.zmid(); if(extrapST.debug() > 2)std::cout << "extrapolating to ST " << std::endl; do { ktrk.extrapolate(tdir,extrapST); @@ -228,16 +230,19 @@ namespace mu2e { } 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){ + double endz = ftraj.position3(endtime).Z() - extrapST.zmid(); // z relative to the ST center + if(endz * startz < 0.0){ + if(extrapST.debug() > 1)std::cout <<"Exiting opposite side ST" << std::endl; return true; + } else { + if(extrapST.debug() > 1)std::cout <<"Exiting same side ST" << std::endl; + return false; } - return false; } template bool KKExtrap::extrapolateTracker(KKTrack& ktrk,TimeDir tdir) const { GeomHandle kkg_h; - ExtrapolateToZ trackerFront(maxdt_,maxdtstep_,btol_,kkg_h->tracker()->front().center().Z(),debug_); + ExtrapolateToZ trackerFront(maxdt_,maxdtstep_,btol_,intertol_,kkg_h->tracker()->front().center().Z(),debug_); if(trackerFront.debug() > 2)std::cout << "extrapolating to Tracker " << std::endl; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId TrackerSID("TT_Front"); @@ -254,7 +259,7 @@ namespace mu2e { template bool KKExtrap::extrapolateTSDA(KKTrack& ktrk,TimeDir tdir) const { GeomHandle kkg_h; - ExtrapolateToZ TSDA(maxdt_,maxdtstep_,btol_,kkg_h->DS()->upstreamAbsorber().center().Z(),debug_); + ExtrapolateToZ TSDA(maxdt_,maxdtstep_,btol_,intertol_,kkg_h->DS()->upstreamAbsorber().center().Z(),debug_); if(TSDA.debug() > 2)std::cout << "extrapolating to TSDA " << std::endl; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId TSDASID("TSDA"); diff --git a/Mu2eKinKal/src/RegrowLoopHelix_module.cc b/Mu2eKinKal/src/RegrowLoopHelix_module.cc index 51fcff5c45..649e016d4c 100644 --- a/Mu2eKinKal/src/RegrowLoopHelix_module.cc +++ b/Mu2eKinKal/src/RegrowLoopHelix_module.cc @@ -9,7 +9,7 @@ #include "fhiclcpp/types/OptionalTable.h" #include "fhiclcpp/types/Tuple.h" #include "fhiclcpp/types/OptionalAtom.h" -#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/EDFilter.h" #include "art/Framework/Principal/Event.h" #include "art/Framework/Principal/Run.h" #include "art/Framework/Principal/Handle.h" @@ -95,13 +95,14 @@ namespace mu2e { fhicl::Table kkfitSettings { Name("KKFitSettings") }; fhicl::Table fitSettings { Name("RefitSettings") }; fhicl::Atom extend {Name("Extend"), Comment("Extend the fit") }; + fhicl::Atom mustRegrow{ Name("MustRegrow"), Comment("Required track regrowing success for event to pass")}; fhicl::OptionalTable extrapSettings { Name("ExtrapolationSettings") }; }; - class RegrowLoopHelix : public art::EDProducer { + class RegrowLoopHelix : public art::EDFilter { public: - using Parameters = art::EDProducer::Table; + using Parameters = art::EDFilter::Table; using KTRAJ = KinKal::LoopHelix; using PKTRAJ = KinKal::ParticleTrajectory; using PKTRAJPTR = std::unique_ptr; @@ -136,9 +137,11 @@ namespace mu2e { using DOMAINPTR = std::shared_ptr; using DOMAINCOL = std::set; + enum regrown {none=0,any,all}; // states for event testing + explicit RegrowLoopHelix(const Parameters& settings); - void beginRun(art::Run& run) override; - void produce(art::Event& event) override; + bool beginRun(art::Run& run) override; + bool filter(art::Event& event) override; void endJob() override; private: int debug_; @@ -158,9 +161,10 @@ namespace mu2e { bool has_cccol_, has_kseedcol_, has_kseedptrcol_; std::unique_ptr extrap_; PDGCode::type tpart_ = PDGCode::unknown; + regrown mustregrow_; // requirement on tracks to regrow for the event to pass }; - RegrowLoopHelix::RegrowLoopHelix(const Parameters& settings) : art::EDProducer(settings), + RegrowLoopHelix::RegrowLoopHelix(const Parameters& settings) : art::EDFilter(settings), debug_(settings().debug()), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), kkfit_(settings().kkfitSettings()), @@ -192,16 +196,26 @@ namespace mu2e { if(has_cccol_) cccol_T_ = art::ProductToken(consumes(settings().caloClusterCollection().value())); int type; if(settings().fitParticle(type))tpart_ = static_cast(type); + + if(settings().mustRegrow() == "None") + mustregrow_ = none; + else if(settings().mustRegrow() == "Any") + mustregrow_ = any; + else if(settings().mustRegrow() == "All") + mustregrow_ = all; + else + throw cet::exception("RECO")<<"mu2e::RegrowLoopHelix: invalid regrow success requirement " << settings().mustRegrow() << " : Enter 'None', 'Any' or 'All'" << endl; } - void RegrowLoopHelix::beginRun(art::Run& run) + bool RegrowLoopHelix::beginRun(art::Run& run) { GeomHandle bfmgr; GeomHandle det; kkbf_ = std::move(std::make_unique(*bfmgr,*det)); + return true; } - void RegrowLoopHelix::produce(art::Event& event) + bool RegrowLoopHelix::filter(art::Event& event) { // proditions auto const& strawresponse = strawResponse_h_.getPtr(event.id()); @@ -249,8 +263,9 @@ namespace mu2e { for (auto const& kseedptr : kseedptrs) { auto const& kseed = *kseedptr; if(!kseed.loopHelixFit())throw cet::exception("RECO")<<"mu2e::RegrowLoopHelix: passed KalSeed from non-LoopHelix fit " << endl; + auto tpart = tpart_ != PDGCode::unknown ? static_cast(std::copysign(static_cast(tpart_),static_cast(kseed.particle()))) : kseed.particle(); // optionally override the particle type // regrow the components from the seed - PKTRAJPTR trajptr = kseed.loopHelixFitTrajectory(tpart_); + PKTRAJPTR trajptr = kseed.loopHelixFitTrajectory(tpart); KKSTRAWHITCOL strawhits; strawhits.reserve(kseed.hits().size()); KKSTRAWXINGCOL strawxings; @@ -276,7 +291,6 @@ namespace mu2e { // require hits and consistent BField domains if(goodhits && (domains.size() > 0 || !config_.bfcorr_)){ // create the KKTrack from these - auto tpart = tpart_ != PDGCode::unknown ? static_cast(std::copysign(static_cast(tpart_),static_cast(kseed.particle()))) : kseed.particle(); // optionally override the particle type auto ktrk = std::make_unique(config_,*kkbf_,tpart,trajptr,strawhits,strawxings,calohits,paramhits, domains); if(ktrk && ktrk->fitStatus().usable()){ if(debug_ > 0) std::cout << "RegrowLoopHelix: successful track refit" << std::endl; @@ -327,11 +341,18 @@ namespace mu2e { } } } + // select events based on regrowing success + bool accept(true); + if(mustregrow_ == all) + accept = ktrkcol->size() == kseedptrs.size(); + else if(mustregrow_ == any) + accept = ktrkcol->size() > 0; // store output event.put(move(ktrkcol)); event.put(move(rgkseedcol)); if(fillMCAssns_)event.put(move(rgksmca)); if(copyKalSeedMCs_)event.put(move(rgksmcc)); + return accept; } void RegrowLoopHelix::endJob() diff --git a/TrkReco/fcl/prolog.fcl b/TrkReco/fcl/prolog.fcl index ed8043143b..8f060b2cbf 100644 --- a/TrkReco/fcl/prolog.fcl +++ b/TrkReco/fcl/prolog.fcl @@ -140,7 +140,7 @@ TrkReco: { @table::TrkReco module_type : SelectReflections Surface : TT_Front MaxDeltaT : 5 # ns - MaxDeltaP : 15 # MeV + MaxDeltaP : 10 # MeV # optionally only merge selected candidates # a simple selector is used here, but any selector implemented as a tool can be used Selector : { @@ -149,7 +149,7 @@ TrkReco: { @table::TrkReco MaxMomentum : 300.0 #MeV/c MinFitConsistency : 1e-6 MinDeltaNHitFraction : 0.0 - MinActiveHits : 0 + MinActiveHits : 15 } SelectBestPair : 3 # select pair with highest # of active hits debugLevel : 0 diff --git a/Validation/src/ValKalSeed.cc b/Validation/src/ValKalSeed.cc index d677813379..56fa0f8f44 100644 --- a/Validation/src/ValKalSeed.cc +++ b/Validation/src/ValKalSeed.cc @@ -164,7 +164,7 @@ namespace mu2e { p = mom3.R(); rho = ikinter->position3().Rho(); } - auto stinters = ks.intersections(SurfaceId("ST_Foils",-1)); // match all foils + auto stinters = ks.intersections(SurfaceId("ST_Foils",SurfaceId::allIndices_)); // match all foils _hNST->Fill(stinters.size()); for (auto stinter : stinters) _hSTdP->Fill(stinter->dMom());