Merge pull request #1595 from brownd1978/kseed

Support for regrowing KKTrack fits from KalSeed objects

Author: brownd1978

.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

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)

DataProducts/inc/IndexMap.hh

@@ -48,6 +48,7 @@ namespace mu2e {
         ost << std::endl;
       }
     }
+    auto const& map() const { return _theMap; }
 
   private:
     std::map<FullIndex, CondensedIndex> _theMap; // FullIndex is first, CondensedIndex is second

EventDisplay/src/DataInterface.cc

@@ -1003,7 +1003,7 @@ void DataInterface::fillEvent(boost::shared_ptr<ContentSelector> const &contentS
       {
         const mu2e::TrkStrawHitSeed &hit = hits.at(j);
         int    sid = hit.strawId().asUint16();
-        double time = hit.hitTime();
+        double time = hit.time();
 
         std::map<int,boost::shared_ptr<Straw> >::iterator straw=_straws.find(sid);
         if(straw!=_straws.end() && !std::isnan(time))

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)

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

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;
           }

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

Mu2eKinKal/inc/KKCaloHit.hh

@@ -26,6 +26,29 @@ namespace mu2e {
       using CA = KinKal::ClosestApproach<KTRAJ,SensorLine>;
       using HIT = KinKal::Hit<KTRAJ>;
       using KTRAJPTR = std::shared_ptr<KTRAJ>;
+      // clone op for reinstantiation
+      KKCaloHit(KKCaloHit<KTRAJ> 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<KTRAJ> > clone(CloneContext& context) const override{
+        auto rv = std::make_shared< KKCaloHit<KTRAJ> >(*this);
+        auto ca = rv->closestApproach();
+        auto trajectory = std::make_shared<KTRAJ>(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 <class KTRAJ> KKCaloHit<KTRAJ>::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 <class KTRAJ> void KKCaloHit<KTRAJ>::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);