Skip to content

Commit dde444e

Browse files
oksuzianoksuzian
authored
Merge pull request #1877 from brownd1978/refreco
Add ability to change particle type when regrowing
2 parents 650ced2 + 23ef8e1 commit dde444e

10 files changed

Lines changed: 89 additions & 31 deletions

File tree

‎Mu2eKinKal/src/RegrowLoopHelix_module.cc‎

Lines changed: 9 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -91,6 +91,7 @@ namespace mu2e {
9191
fhicl::Atom<art::InputTag> indexMap {Name("StrawDigiIndexMap"), Comment("Map between original and reduced ComboHits") };
9292
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") };
9393
fhicl::OptionalAtom<bool> copyKalSeedMCs { Name("CopyKalSeedMCs"), Comment("If set, and KalSeedMCs are referenced, create a deep copy of the input") };
94+
fhicl::OptionalAtom<int> fitParticle { Name("FitParticle"), Comment("PDG code of particle type to refit: e-, e+, mu-, ...")};
9495
fhicl::Table<KKFitConfig> kkfitSettings { Name("KKFitSettings") };
9596
fhicl::Table<KKConfig> fitSettings { Name("RefitSettings") };
9697
fhicl::Atom<bool> extend {Name("Extend"), Comment("Extend the fit") };
@@ -152,10 +153,11 @@ namespace mu2e {
152153
art::ProductToken<CaloClusterCollection> cccol_T_;
153154
art::ProductToken<IndexMap> indexmap_T_;
154155
art::InputTag ksmca_T_;
155-
bool fillMCAssns_, copyKalSeedMCs_;
156+
bool fillMCAssns_ = false, copyKalSeedMCs_ = false;
156157
bool extend_;
157158
bool has_cccol_, has_kseedcol_, has_kseedptrcol_;
158159
std::unique_ptr<KKExtrap> extrap_;
160+
PDGCode::type tpart_ = PDGCode::unknown;
159161
};
160162

161163
RegrowLoopHelix::RegrowLoopHelix(const Parameters& settings) : art::EDProducer(settings),
@@ -188,6 +190,8 @@ namespace mu2e {
188190
if(has_kseedcol_)kseedcol_T_ = art::ProductToken<KalSeedCollection>(consumes<KalSeedCollection>(settings().kalSeedCollection().value()));
189191
if(has_kseedptrcol_)kseedptrcol_T_ = art::ProductToken<KalSeedPtrCollection>(consumes<KalSeedPtrCollection>(settings().kalSeedPtrCollection().value()));
190192
if(has_cccol_) cccol_T_ = art::ProductToken<CaloClusterCollection>(consumes<CaloClusterCollection>(settings().caloClusterCollection().value()));
193+
int type;
194+
if(settings().fitParticle(type))tpart_ = static_cast<PDGCode::type>(type);
191195
}
192196

193197
void RegrowLoopHelix::beginRun(art::Run& run)
@@ -246,7 +250,7 @@ namespace mu2e {
246250
auto const& kseed = *kseedptr;
247251
if(!kseed.loopHelixFit())throw cet::exception("RECO")<<"mu2e::RegrowLoopHelix: passed KalSeed from non-LoopHelix fit " << endl;
248252
// regrow the components from the seed
249-
PKTRAJPTR trajptr = kseed.loopHelixFitTrajectory();
253+
PKTRAJPTR trajptr = kseed.loopHelixFitTrajectory(tpart_);
250254
KKSTRAWHITCOL strawhits;
251255
strawhits.reserve(kseed.hits().size());
252256
KKSTRAWXINGCOL strawxings;
@@ -271,8 +275,9 @@ namespace mu2e {
271275
if(debug_ > 5)static_cast<KinKal::PiecewiseTrajectory<KTRAJ>*>(trajptr.get())->print(std::cout,2);
272276
// require hits and consistent BField domains
273277
if(goodhits && (domains.size() > 0 || !config_.bfcorr_)){
274-
// create the KKTrack from these components
275-
auto ktrk = std::make_unique<KKTRK>(config_,*kkbf_,kseed.particle(),trajptr,strawhits,strawxings,calohits,paramhits,domains);
278+
// create the KKTrack from these
279+
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
280+
auto ktrk = std::make_unique<KKTRK>(config_,*kkbf_,tpart,trajptr,strawhits,strawxings,calohits,paramhits, domains);
276281
if(ktrk && ktrk->fitStatus().usable()){
277282
if(debug_ > 0) std::cout << "RegrowLoopHelix: successful track refit" << std::endl;
278283
if(extend_)kkfit_.extendTrack(config_,*kkbf_, *tracker,*strawresponse, chcol, *calo_h, cc_H , *ktrk );

‎Print/fcl/printPrimaryParticle.fcl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -18,14 +18,14 @@ physics :{
1818
printModule : {
1919
module_type : PrintModule
2020
primaryParticlePrinter : {
21-
verbose : 2
21+
verbose : 1
2222
}
2323
} # printModule
2424

2525

2626
} # analyzers
2727

28-
ana : [ printModule, printProductList ]
28+
ana : [ printModule]
2929
end_paths : [ ana ]
3030

3131
}

‎Print/inc/PrimaryParticlePrinter.hh‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,10 @@ namespace mu2e {
1919
class PrimaryParticlePrinter : public ProductPrinter {
2020
public:
2121
PrimaryParticlePrinter() {}
22-
PrimaryParticlePrinter(const Config& conf) : ProductPrinter(conf) {}
22+
PrimaryParticlePrinter(const Config& conf) : ProductPrinter(conf) {
23+
gprint_.setVerbose(verbose());
24+
sprint_.setVerbose(verbose());
25+
}
2326

2427
// all the ways to request a printout
2528
void Print(art::Event const& event, std::ostream& os = std::cout) override;

‎RecoDataProducts/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,7 @@ cet_make_library(
3535

3636
Offline::DataProducts
3737
Offline::GeneralUtilities
38+
Offline::GlobalConstantsService
3839
Offline::KinKalGeom
3940
Offline::Mu2eKinKal
4041
Offline::TrackerConditions

‎RecoDataProducts/inc/KalSeed.hh‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -71,9 +71,10 @@ namespace mu2e {
7171
// reconstitute (as best as possible) the fit trajectory. The ptr will be null if the fit wasn't based on the requested trajector type
7272
// Note these return by value
7373
// Note that the returned piecetraj may have large gaps, unless the full fit trajectory was stored in the seed.
74-
LHPTPtr loopHelixFitTrajectory() const;
75-
CHPTPtr centralHelixFitTrajectory() const;
76-
KLPTPtr kinematicLineFitTrajectory() const;
74+
// Optionally override the particle Id
75+
LHPTPtr loopHelixFitTrajectory(PDGCode::type ptype = PDGCode::unknown) const;
76+
CHPTPtr centralHelixFitTrajectory(PDGCode::type ptype = PDGCode::unknown) const;
77+
KLPTPtr kinematicLineFitTrajectory(PDGCode::type ptype = PDGCode::unknown) const;
7778

7879
// global information about the track
7980
PDGCode::type _tpart = PDGCode::unknown; // particle assumed for this fit

‎RecoDataProducts/inc/KalSegment.hh‎

Lines changed: 4 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -35,12 +35,10 @@ namespace mu2e {
3535
KinKal::VEC3 momentum3() const { return _pstate.momentum3(); }
3636
KinKal::VEC3 velocity() const { return _pstate.velocity(); }
3737
KinKal::VEC3 position3() const { return _pstate.position3(); }
38-
// convert content to a LoopHelix
39-
KinKal::LoopHelix loopHelix() const { return KinKal::LoopHelix(_pstate, bnom(),timeRange()); }
40-
// convert to a CentralHelix
41-
KinKal::CentralHelix centralHelix() const { return KinKal::CentralHelix(_pstate, bnom(),timeRange()); }
42-
// convert to a KinematicLine
43-
KinKal::KinematicLine kinematicLine() const { return KinKal::KinematicLine(_pstate, bnom(),timeRange()); }
38+
// convert content to a Trajectory. Optionally override the particle mass
39+
KinKal::LoopHelix loopHelix(double mass = -1.0) const;
40+
KinKal::CentralHelix centralHelix(double mass = -1.0) const;
41+
KinKal::KinematicLine kinematicLine(double mass = -1.0) const;
4442
double const& tmin() const { return _tmin; }
4543
double const& tmax() const { return _tmax; }
4644
KinKal::TimeRange timeRange() const;

‎RecoDataProducts/src/KalSeed.cc‎

Lines changed: 30 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,26 @@
11
#include "Offline/RecoDataProducts/inc/KalSeed.hh"
2+
#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh"
3+
#include "Offline/GlobalConstantsService/inc/ParticleDataList.hh"
4+
#include "cetlib_except/exception.h"
25
#include <limits>
36
namespace mu2e {
47

58
const double KalSeed::_regrowtol(1e-3); // 1 ps minimum for regrown segments
69

7-
KalSeed::LHPTPtr KalSeed::loopHelixFitTrajectory() const {
10+
KalSeed::LHPTPtr KalSeed::loopHelixFitTrajectory(PDGCode::type ptype) const {
11+
double mass = -1.0;
12+
if(ptype != PDGCode::unknown && ptype != _tpart){
13+
auto const& ptable = GlobalConstantsHandle<ParticleDataList>();
14+
mass = ptable->particle(ptype).mass();
15+
}
16+
817
if(loopHelixFit() && segments().size() > 0){
918
// initialize the piecewise trajectory with the front segment
10-
LHPTPtr ptraj(new KalSeed::LHPT(segments().front().loopHelix()));
19+
LHPTPtr ptraj(new KalSeed::LHPT(segments().front().loopHelix(mass)));
1120
auto iseg = segments().begin(); ++iseg;
1221
while(iseg != segments().end()){
1322
if(iseg->timeRange().range() > _regrowtol ){
14-
auto trajptr = std::make_shared<KinKal::LoopHelix>(iseg->loopHelix());
23+
auto trajptr = std::make_shared<KinKal::LoopHelix>(iseg->loopHelix(mass));
1524
ptraj->add(trajptr); // note this call resolves the phi0 ambiguity
1625
}
1726
++iseg;
@@ -21,14 +30,20 @@ namespace mu2e {
2130
return LHPTPtr();
2231
}
2332

24-
KalSeed::CHPTPtr KalSeed::centralHelixFitTrajectory() const {
33+
KalSeed::CHPTPtr KalSeed::centralHelixFitTrajectory(PDGCode::type ptype) const {
34+
double mass = -1.0;
35+
if(ptype != PDGCode::unknown && ptype != _tpart){
36+
auto const& ptable = GlobalConstantsHandle<ParticleDataList>();
37+
mass = ptable->particle(ptype).mass();
38+
}
39+
2540
if(centralHelixFit() && segments().size() > 0){
2641
// initialize the piecewise trajectory with the front segment
27-
CHPTPtr ptraj(new KalSeed::CHPT(segments().front().centralHelix()));
42+
CHPTPtr ptraj(new KalSeed::CHPT(segments().front().centralHelix(mass)));
2843
auto iseg = segments().begin(); ++iseg;
2944
while(iseg != segments().end()){
3045
if(iseg->timeRange().range() > _regrowtol ){
31-
auto trajptr = std::make_shared<KinKal::CentralHelix>(iseg->centralHelix());
46+
auto trajptr = std::make_shared<KinKal::CentralHelix>(iseg->centralHelix(mass));
3247
ptraj->add(trajptr);
3348
}
3449
++iseg;
@@ -38,14 +53,20 @@ namespace mu2e {
3853
return CHPTPtr();
3954
}
4055

41-
KalSeed::KLPTPtr KalSeed::kinematicLineFitTrajectory() const {
56+
KalSeed::KLPTPtr KalSeed::kinematicLineFitTrajectory(PDGCode::type ptype) const {
57+
double mass = -1.0;
58+
if(ptype != PDGCode::unknown && ptype != _tpart){
59+
auto const& ptable = GlobalConstantsHandle<ParticleDataList>();
60+
mass = ptable->particle(ptype).mass();
61+
}
62+
4263
if(kinematicLineFit() && segments().size() > 0){
4364
// initialize the piecewise trajectory with the front segment
44-
KLPTPtr ptraj(new KalSeed::KLPT(segments().front().kinematicLine()));
65+
KLPTPtr ptraj(new KalSeed::KLPT(segments().front().kinematicLine(mass)));
4566
auto iseg = segments().begin(); ++iseg;
4667
while(iseg != segments().end()){
4768
if(iseg->timeRange().range() > _regrowtol ){
48-
auto trajptr = std::make_shared<KinKal::KinematicLine>(iseg->kinematicLine());
69+
auto trajptr = std::make_shared<KinKal::KinematicLine>(iseg->kinematicLine(mass));
4970
ptraj->add(trajptr);
5071
}
5172
++iseg;

‎RecoDataProducts/src/KalSegment.cc‎

Lines changed: 32 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,38 @@
11
#include "Offline/RecoDataProducts/inc/KalSegment.hh"
22
#include "cetlib_except/exception.h"
3-
43
namespace mu2e {
54

5+
KinKal::LoopHelix KalSegment::loopHelix(double mass) const {
6+
if(mass < 0.0){
7+
return KinKal::LoopHelix(_pstate, bnom(),timeRange());
8+
} else {
9+
return KinKal::LoopHelix(KinKal::ParticleStateEstimate(
10+
KinKal::ParticleState(_pstate.position3(),_pstate.momentum3(),_pstate.time(),mass,_pstate.charge()),
11+
_pstate.stateCovariance()),
12+
bnom(),timeRange());
13+
}
14+
}
15+
KinKal::CentralHelix KalSegment::centralHelix(double mass ) const {
16+
if(mass < 0.0){
17+
return KinKal::CentralHelix(_pstate, bnom(),timeRange());
18+
} else {
19+
return KinKal::CentralHelix(KinKal::ParticleStateEstimate(
20+
KinKal::ParticleState(_pstate.position3(),_pstate.momentum3(),_pstate.time(),mass,_pstate.charge()),
21+
_pstate.stateCovariance()),
22+
bnom(),timeRange());
23+
}
24+
}
25+
KinKal::KinematicLine KalSegment::kinematicLine(double mass ) const {
26+
if(mass < 0.0){
27+
return KinKal::KinematicLine(_pstate, bnom(),timeRange());
28+
} else {
29+
return KinKal::KinematicLine(KinKal::ParticleStateEstimate(
30+
KinKal::ParticleState(_pstate.position3(),_pstate.momentum3(),_pstate.time(),mass,_pstate.charge()),
31+
_pstate.stateCovariance()),
32+
bnom(),timeRange());
33+
}
34+
}
35+
636
KinKal::TimeRange KalSegment::timeRange() const {
737
// protect against unphysical times
838
if(_tmin <= _tmax){
@@ -25,10 +55,7 @@ namespace mu2e {
2555
auto kl = kinematicLine();
2656
return kl.paramVal(KinKal::KinematicLine::t0Index());
2757
}
28-
// throw cet::exception("RECO")<<"mu2e::KalSegment: no trajectory specified in flag" << std::endl;
29-
// for now, revert to the legacy implementation. Once BTrk is fully removed this should be removed
30-
auto vel = _pstate.velocity();
31-
return _pstate.time() - _pstate.position3().Z()/vel.Z();
58+
throw cet::exception("RECO")<<"mu2e::KalSegment: no trajectory specified in flag" << std::endl;
3259
}
3360

3461
double KalSegment::t0Var(TrkFitFlag const& flag) const {

‎RecoDataProducts/src/SConscript‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,8 +16,10 @@ mainlib = helper.make_mainlib ( [ 'mu2e_DataProducts',
1616
'mu2e_GeomPrimitives',
1717
'mu2e_GeometryService',
1818
'mu2e_GeneralUtilities',
19+
'mu2e_GlobalConstantsService',
1920
'art_Framework_Core',
2021
'art_Framework_Principal',
22+
'art_Framework_Services_Registry',
2123
'art_Persistency_Common',
2224
'art_Persistency_Provenance',
2325
'art_Utilities',

‎TrkReco/fcl/prolog.fcl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -139,7 +139,7 @@ TrkReco: { @table::TrkReco
139139
SelectReflections: {
140140
module_type : SelectReflections
141141
Surface : TT_Front
142-
MaxDeltaT : 15 # ns
142+
MaxDeltaT : 5 # ns
143143
MaxDeltaP : 15 # MeV
144144
# optionally only merge selected candidates
145145
# a simple selector is used here, but any selector implemented as a tool can be used

0 commit comments

Comments
 (0)