diff --git a/Mu2eKinKal/src/RegrowLoopHelix_module.cc b/Mu2eKinKal/src/RegrowLoopHelix_module.cc index f6d6aabfb3..51fcff5c45 100644 --- a/Mu2eKinKal/src/RegrowLoopHelix_module.cc +++ b/Mu2eKinKal/src/RegrowLoopHelix_module.cc @@ -91,6 +91,7 @@ namespace mu2e { 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::OptionalAtom copyKalSeedMCs { Name("CopyKalSeedMCs"), Comment("If set, and KalSeedMCs are referenced, create a deep copy of the input") }; + fhicl::OptionalAtom fitParticle { Name("FitParticle"), Comment("PDG code of particle type to refit: e-, e+, mu-, ...")}; fhicl::Table kkfitSettings { Name("KKFitSettings") }; fhicl::Table fitSettings { Name("RefitSettings") }; fhicl::Atom extend {Name("Extend"), Comment("Extend the fit") }; @@ -152,10 +153,11 @@ namespace mu2e { art::ProductToken cccol_T_; art::ProductToken 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 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(consumes(settings().kalSeedCollection().value())); if(has_kseedptrcol_)kseedptrcol_T_ = art::ProductToken(consumes(settings().kalSeedPtrCollection().value())); if(has_cccol_) cccol_T_ = art::ProductToken(consumes(settings().caloClusterCollection().value())); + int type; + if(settings().fitParticle(type))tpart_ = static_cast(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*>(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,paramhits,domains); + // 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; if(extend_)kkfit_.extendTrack(config_,*kkbf_, *tracker,*strawresponse, chcol, *calo_h, cc_H , *ktrk ); diff --git a/Print/fcl/printPrimaryParticle.fcl b/Print/fcl/printPrimaryParticle.fcl index ded50ef9ed..0e2b055e3b 100644 --- a/Print/fcl/printPrimaryParticle.fcl +++ b/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 ] } diff --git a/Print/inc/PrimaryParticlePrinter.hh b/Print/inc/PrimaryParticlePrinter.hh index ee77c663f3..5f13d5ec43 100644 --- a/Print/inc/PrimaryParticlePrinter.hh +++ b/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; diff --git a/RecoDataProducts/CMakeLists.txt b/RecoDataProducts/CMakeLists.txt index 1a718875fe..6ddb9ac3b0 100644 --- a/RecoDataProducts/CMakeLists.txt +++ b/RecoDataProducts/CMakeLists.txt @@ -35,6 +35,7 @@ cet_make_library( Offline::DataProducts Offline::GeneralUtilities + Offline::GlobalConstantsService Offline::KinKalGeom Offline::Mu2eKinKal Offline::TrackerConditions diff --git a/RecoDataProducts/inc/KalSeed.hh b/RecoDataProducts/inc/KalSeed.hh index ab01e5db8f..f229e33cb6 100644 --- a/RecoDataProducts/inc/KalSeed.hh +++ b/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 diff --git a/RecoDataProducts/inc/KalSegment.hh b/RecoDataProducts/inc/KalSegment.hh index da48160d3f..593c218e8b 100644 --- a/RecoDataProducts/inc/KalSegment.hh +++ b/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; diff --git a/RecoDataProducts/src/KalSeed.cc b/RecoDataProducts/src/KalSeed.cc index 16a8432e73..a2314c400d 100644 --- a/RecoDataProducts/src/KalSeed.cc +++ b/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 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(); + 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(iseg->loopHelix()); + auto trajptr = std::make_shared(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(); + 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(iseg->centralHelix()); + auto trajptr = std::make_shared(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(); + 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(iseg->kinematicLine()); + auto trajptr = std::make_shared(iseg->kinematicLine(mass)); ptraj->add(trajptr); } ++iseg; diff --git a/RecoDataProducts/src/KalSegment.cc b/RecoDataProducts/src/KalSegment.cc index 804e67ad2f..bb908f192c 100644 --- a/RecoDataProducts/src/KalSegment.cc +++ b/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 { diff --git a/RecoDataProducts/src/SConscript b/RecoDataProducts/src/SConscript index 62c9fa3b97..65a2f8a2f0 100644 --- a/RecoDataProducts/src/SConscript +++ b/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', diff --git a/TrkReco/fcl/prolog.fcl b/TrkReco/fcl/prolog.fcl index 6f6f1f8c64..ed8043143b 100644 --- a/TrkReco/fcl/prolog.fcl +++ b/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