Add ability to change particle type when regrowing

Mu2eKinKal/src/RegrowLoopHelix_module.cc

@@ -91,6 +91,7 @@ namespace mu2e {
     fhicl::Atom<art::InputTag> indexMap {Name("StrawDigiIndexMap"), Comment("Map between original and reduced ComboHits") };
     fhicl::OptionalAtom<art::InputTag> kalSeedMCAssns {Name("KalSeedMCAssns"), Comment("Association to KalSeedMC. If set, regrown KalSeeds will be associated with the same KalSeedMC as the original") };
     fhicl::OptionalAtom<bool> copyKalSeedMCs { Name("CopyKalSeedMCs"), Comment("If set, and KalSeedMCs are referenced, create a deep copy of the input") };
+    fhicl::OptionalAtom<int> fitParticle {  Name("FitParticle"), Comment("PDG code of particle type to refit: e-, e+, mu-, ...")};
     fhicl::Table<KKFitConfig> kkfitSettings { Name("KKFitSettings") };
     fhicl::Table<KKConfig> fitSettings { Name("RefitSettings") };
     fhicl::Atom<bool> extend {Name("Extend"), Comment("Extend the fit") };
@@ -152,10 +153,11 @@ namespace mu2e {
       art::ProductToken<CaloClusterCollection> cccol_T_;
       art::ProductToken<IndexMap> indexmap_T_;
       art::InputTag ksmca_T_;
-      bool fillMCAssns_, copyKalSeedMCs_;
+      bool fillMCAssns_ = false, copyKalSeedMCs_ = false;
       bool extend_;
       bool has_cccol_, has_kseedcol_, has_kseedptrcol_;
       std::unique_ptr<KKExtrap> extrap_;
+      PDGCode::type tpart_ = PDGCode::unknown;
   };
 
   RegrowLoopHelix::RegrowLoopHelix(const Parameters& settings) : art::EDProducer(settings),
@@ -188,6 +190,8 @@ namespace mu2e {
     if(has_kseedcol_)kseedcol_T_ = art::ProductToken<KalSeedCollection>(consumes<KalSeedCollection>(settings().kalSeedCollection().value()));
     if(has_kseedptrcol_)kseedptrcol_T_ = art::ProductToken<KalSeedPtrCollection>(consumes<KalSeedPtrCollection>(settings().kalSeedPtrCollection().value()));
     if(has_cccol_) cccol_T_ = art::ProductToken<CaloClusterCollection>(consumes<CaloClusterCollection>(settings().caloClusterCollection().value()));
+    int type;
+    if(settings().fitParticle(type))tpart_ = static_cast<PDGCode::type>(type);
   }
 
   void RegrowLoopHelix::beginRun(art::Run& run)
@@ -246,7 +250,7 @@ namespace mu2e {
       auto const& kseed = *kseedptr;
       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();
+      PKTRAJPTR trajptr = kseed.loopHelixFitTrajectory(tpart_);
       KKSTRAWHITCOL strawhits;
       strawhits.reserve(kseed.hits().size());
       KKSTRAWXINGCOL strawxings;
@@ -271,8 +275,9 @@ namespace mu2e {
       if(debug_ > 5)static_cast<KinKal::PiecewiseTrajectory<KTRAJ>*>(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<KKTRK>(config_,*kkbf_,kseed.particle(),trajptr,strawhits,strawxings,calohits,paramhits,domains);
+      // create the KKTrack from these
+        auto tpart = tpart_ != PDGCode::unknown ? static_cast<PDGCode>(std::copysign(static_cast<int>(tpart_),static_cast<int>(kseed.particle()))) : kseed.particle(); // optionally override the particle type
+        auto ktrk = std::make_unique<KKTRK>(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;
           if(extend_)kkfit_.extendTrack(config_,*kkbf_, *tracker,*strawresponse, chcol, *calo_h, cc_H , *ktrk );

Print/fcl/printPrimaryParticle.fcl

@@ -18,14 +18,14 @@ physics :{
     printModule : {
       module_type : PrintModule
       primaryParticlePrinter : {
-        verbose : 2
+        verbose : 1
       }
     } # printModule
 
 
   }  # analyzers
 
-  ana       : [ printModule, printProductList ]
+  ana       : [ printModule]
   end_paths : [ ana ]
 
 }

Print/inc/PrimaryParticlePrinter.hh

@@ -19,7 +19,10 @@ namespace mu2e {
 class PrimaryParticlePrinter : public ProductPrinter {
  public:
   PrimaryParticlePrinter() {}
-  PrimaryParticlePrinter(const Config& conf) : ProductPrinter(conf) {}
+  PrimaryParticlePrinter(const Config& conf) : ProductPrinter(conf) {
+    gprint_.setVerbose(verbose());
+    sprint_.setVerbose(verbose());
+  }
 
   // all the ways to request a printout
   void Print(art::Event const& event, std::ostream& os = std::cout) override;

RecoDataProducts/CMakeLists.txt

@@ -35,6 +35,7 @@ cet_make_library(
 
       Offline::DataProducts
       Offline::GeneralUtilities
+      Offline::GlobalConstantsService
       Offline::KinKalGeom
       Offline::Mu2eKinKal
       Offline::TrackerConditions

RecoDataProducts/inc/KalSeed.hh

@@ -71,9 +71,10 @@ namespace mu2e {
     // reconstitute (as best as possible) the fit trajectory.  The ptr will be null if the fit wasn't based on the requested trajector type
     // Note these return by value
     // Note that the returned piecetraj may have large gaps, unless the full fit trajectory was stored in the seed.
-    LHPTPtr loopHelixFitTrajectory() const;
-    CHPTPtr centralHelixFitTrajectory() const;
-    KLPTPtr kinematicLineFitTrajectory() const;
+    // Optionally override the particle Id
+    LHPTPtr loopHelixFitTrajectory(PDGCode::type ptype = PDGCode::unknown) const;
+    CHPTPtr centralHelixFitTrajectory(PDGCode::type ptype = PDGCode::unknown) const;
+    KLPTPtr kinematicLineFitTrajectory(PDGCode::type ptype = PDGCode::unknown) const;
 
     // global information about the track
     PDGCode::type   _tpart = PDGCode::unknown; // particle assumed for this fit

RecoDataProducts/inc/KalSegment.hh

@@ -35,12 +35,10 @@ namespace mu2e {
     KinKal::VEC3 momentum3() const { return _pstate.momentum3(); }
     KinKal::VEC3 velocity() const { return _pstate.velocity(); }
     KinKal::VEC3 position3() const { return _pstate.position3(); }
-    // convert content to a LoopHelix
-    KinKal::LoopHelix loopHelix() const { return KinKal::LoopHelix(_pstate, bnom(),timeRange()); }
-    // convert to a CentralHelix
-    KinKal::CentralHelix centralHelix() const { return KinKal::CentralHelix(_pstate, bnom(),timeRange()); }
-    // convert to a KinematicLine
-    KinKal::KinematicLine kinematicLine() const { return KinKal::KinematicLine(_pstate, bnom(),timeRange()); }
+    // convert content to a Trajectory. Optionally override the particle mass
+    KinKal::LoopHelix loopHelix(double mass = -1.0) const;
+    KinKal::CentralHelix centralHelix(double mass = -1.0) const;
+    KinKal::KinematicLine kinematicLine(double mass = -1.0) const;
     double const& tmin() const { return _tmin; }
     double const& tmax() const { return _tmax; }
     KinKal::TimeRange timeRange() const;

RecoDataProducts/src/KalSeed.cc

@@ -1,17 +1,26 @@
 #include "Offline/RecoDataProducts/inc/KalSeed.hh"
+#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh"
+#include "Offline/GlobalConstantsService/inc/ParticleDataList.hh"
+#include "cetlib_except/exception.h"
 #include <limits>
 namespace mu2e {
 
   const double KalSeed::_regrowtol(1e-3); // 1 ps minimum for regrown segments
 
-  KalSeed::LHPTPtr KalSeed::loopHelixFitTrajectory() const {
+  KalSeed::LHPTPtr KalSeed::loopHelixFitTrajectory(PDGCode::type ptype) const {
+    double mass = -1.0;
+    if(ptype != PDGCode::unknown && ptype != _tpart){
+      auto const& ptable = GlobalConstantsHandle<ParticleDataList>();
+      mass = ptable->particle(ptype).mass();
+    }
+
     if(loopHelixFit() && segments().size() > 0){
       // initialize the piecewise trajectory with the front segment
-      LHPTPtr ptraj(new KalSeed::LHPT(segments().front().loopHelix()));
+      LHPTPtr ptraj(new KalSeed::LHPT(segments().front().loopHelix(mass)));
       auto iseg = segments().begin(); ++iseg;
       while(iseg != segments().end()){
         if(iseg->timeRange().range() > _regrowtol ){
-          auto trajptr = std::make_shared<KinKal::LoopHelix>(iseg->loopHelix());
+          auto trajptr = std::make_shared<KinKal::LoopHelix>(iseg->loopHelix(mass));
           ptraj->add(trajptr); // note this call resolves the phi0 ambiguity
         }
         ++iseg;
@@ -21,14 +30,20 @@ namespace mu2e {
       return LHPTPtr();
   }
 
-  KalSeed::CHPTPtr KalSeed::centralHelixFitTrajectory() const {
+  KalSeed::CHPTPtr KalSeed::centralHelixFitTrajectory(PDGCode::type ptype) const {
+    double mass = -1.0;
+    if(ptype != PDGCode::unknown && ptype != _tpart){
+      auto const& ptable = GlobalConstantsHandle<ParticleDataList>();
+      mass = ptable->particle(ptype).mass();
+    }
+
     if(centralHelixFit() && segments().size() > 0){
       // initialize the piecewise trajectory with the front segment
-      CHPTPtr ptraj(new KalSeed::CHPT(segments().front().centralHelix()));
+      CHPTPtr ptraj(new KalSeed::CHPT(segments().front().centralHelix(mass)));
       auto iseg = segments().begin(); ++iseg;
       while(iseg != segments().end()){
         if(iseg->timeRange().range() > _regrowtol ){
-          auto trajptr = std::make_shared<KinKal::CentralHelix>(iseg->centralHelix());
+          auto trajptr = std::make_shared<KinKal::CentralHelix>(iseg->centralHelix(mass));
           ptraj->add(trajptr);
         }
         ++iseg;
@@ -38,14 +53,20 @@ namespace mu2e {
       return CHPTPtr();
   }
 
-  KalSeed::KLPTPtr KalSeed::kinematicLineFitTrajectory() const {
+  KalSeed::KLPTPtr KalSeed::kinematicLineFitTrajectory(PDGCode::type ptype) const {
+    double mass = -1.0;
+    if(ptype != PDGCode::unknown && ptype != _tpart){
+      auto const& ptable = GlobalConstantsHandle<ParticleDataList>();
+      mass = ptable->particle(ptype).mass();
+    }
+
     if(kinematicLineFit() && segments().size() > 0){
       // initialize the piecewise trajectory with the front segment
-      KLPTPtr ptraj(new KalSeed::KLPT(segments().front().kinematicLine()));
+      KLPTPtr ptraj(new KalSeed::KLPT(segments().front().kinematicLine(mass)));
       auto iseg = segments().begin(); ++iseg;
       while(iseg != segments().end()){
         if(iseg->timeRange().range() > _regrowtol ){
-          auto trajptr = std::make_shared<KinKal::KinematicLine>(iseg->kinematicLine());
+          auto trajptr = std::make_shared<KinKal::KinematicLine>(iseg->kinematicLine(mass));
           ptraj->add(trajptr);
         }
         ++iseg;

RecoDataProducts/src/KalSegment.cc

@@ -1,8 +1,38 @@
 #include "Offline/RecoDataProducts/inc/KalSegment.hh"
 #include "cetlib_except/exception.h"
-
 namespace mu2e {
 
+  KinKal::LoopHelix KalSegment::loopHelix(double mass) const {
+    if(mass < 0.0){
+      return KinKal::LoopHelix(_pstate, bnom(),timeRange());
+    } else {
+      return KinKal::LoopHelix(KinKal::ParticleStateEstimate(
+            KinKal::ParticleState(_pstate.position3(),_pstate.momentum3(),_pstate.time(),mass,_pstate.charge()),
+            _pstate.stateCovariance()),
+          bnom(),timeRange());
+    }
+  }
+  KinKal::CentralHelix KalSegment::centralHelix(double mass ) const {
+    if(mass < 0.0){
+      return KinKal::CentralHelix(_pstate, bnom(),timeRange());
+    } else {
+      return KinKal::CentralHelix(KinKal::ParticleStateEstimate(
+            KinKal::ParticleState(_pstate.position3(),_pstate.momentum3(),_pstate.time(),mass,_pstate.charge()),
+            _pstate.stateCovariance()),
+          bnom(),timeRange());
+    }
+  }
+  KinKal::KinematicLine KalSegment::kinematicLine(double mass ) const {
+    if(mass < 0.0){
+      return KinKal::KinematicLine(_pstate, bnom(),timeRange());
+    } else {
+      return KinKal::KinematicLine(KinKal::ParticleStateEstimate(
+            KinKal::ParticleState(_pstate.position3(),_pstate.momentum3(),_pstate.time(),mass,_pstate.charge()),
+            _pstate.stateCovariance()),
+          bnom(),timeRange());
+    }
+  }
+
   KinKal::TimeRange KalSegment::timeRange() const {
     // protect against unphysical times
     if(_tmin <= _tmax){
@@ -25,10 +55,7 @@ namespace mu2e {
       auto kl = kinematicLine();
       return kl.paramVal(KinKal::KinematicLine::t0Index());
     }
-    //      throw cet::exception("RECO")<<"mu2e::KalSegment: no trajectory specified in flag" << std::endl;
-    // for now, revert to the legacy implementation.  Once BTrk is fully removed this should be removed
-    auto vel = _pstate.velocity();
-    return _pstate.time() - _pstate.position3().Z()/vel.Z();
+    throw cet::exception("RECO")<<"mu2e::KalSegment: no trajectory specified in flag" << std::endl;
   }
 
   double KalSegment::t0Var(TrkFitFlag const& flag) const {

RecoDataProducts/src/SConscript

@@ -16,8 +16,10 @@ mainlib = helper.make_mainlib ( [ 'mu2e_DataProducts',
   'mu2e_GeomPrimitives',
   'mu2e_GeometryService',
   'mu2e_GeneralUtilities',
+  'mu2e_GlobalConstantsService',
   'art_Framework_Core',
   'art_Framework_Principal',
+  'art_Framework_Services_Registry',
   'art_Persistency_Common',
   'art_Persistency_Provenance',
   'art_Utilities',

TrkReco/fcl/prolog.fcl

@@ -139,7 +139,7 @@ TrkReco: { @table::TrkReco
   SelectReflections: {
     module_type : SelectReflections
     Surface : TT_Front
-    MaxDeltaT : 15 # ns
+    MaxDeltaT : 5 # ns
     MaxDeltaP : 15 # MeV
     # optionally only merge selected candidates
     # a simple selector is used here, but any selector implemented as a tool can be used