Improvements for reflection tracking and energy loss

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<enum_type,std::string> 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<enum_type,std::string> const& names();
   };
   using SurfaceIdEnum = EnumToStringSparse<SurfaceIdDetail>;
   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
   };

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<SurfaceIdEnum::enum_type, char const*>;

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

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

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',

Mu2eKinKal/fcl/prolog.fcl

@@ -636,9 +636,10 @@ Mu2eKinKal : {
     }
     ExtrapolationSettings :{
       @table::Mu2eKinKal.LHDRIFTXTRAP
-      BackToTracker: true
+      BackToTracker : true
     }
     Extend: false
+    MustRegrow : "Any"
   }
 
   producers : {

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 <class KTRAJ> bool needsExtrapolation(KinKal::ParticleTrajectory<KTRAJ> const& fittraj, TimeDir tdir) const;
@@ -50,47 +51,41 @@ namespace mu2e {
   };
 
   template <class KTRAJ> bool ExtrapolateIPA::needsExtrapolation(KinKal::ParticleTrajectory<KTRAJ> 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

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 <class KTRAJ> bool needsExtrapolation(KinKal::ParticleTrajectory<KTRAJ> 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<size_t>& 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<size_t> 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<size_t>& 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<size_t>(ifoil);
+        auto const& foilptr = foils_[jfoil];
+        if(foilptr->center().Z() > zmin && foilptr->center().Z() < zmax ){
+          foils.push_back(jfoil);
         }
       }
     }
-    return retval;
   }
 
 }