Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

MTD reconstruction: propagation of track time of flight uncertainty to vertex reconstruction #44562

Merged
merged 6 commits into from
Apr 12, 2024
Merged

MTD reconstruction: propagation of track time of flight uncertainty to vertex reconstruction #44562

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
52f2b56
Propagated track sigma(TOF) on MTD hit-to-track matching, final track…
noepalm Mar 20, 2024
c74aa70
Updated array from C-style array to std::array
noepalm Mar 26, 2024
d69c2bd
Added kMixd calculation strategy for sigma(TOF), set as default in ca…
noepalm Mar 26, 2024
3d658f8
Code checks && format
noepalm Mar 26, 2024
6d4465f
Added dterror2 to track TOF PID info struct and used it to replace oc…
noepalm Apr 8, 2024
d293345
Code checks and format
noepalm Apr 9, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 32 additions & 9 deletions RecoMTD/TimingIDTools/plugins/TOFPIDProducer.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,9 @@ class TOFPIDProducer : public edm::stream::EDProducer<> {
edm::EDGetTokenT<edm::ValueMap<float>> sigmatmtdToken_;
edm::EDGetTokenT<edm::ValueMap<float>> tofkToken_;
edm::EDGetTokenT<edm::ValueMap<float>> tofpToken_;
edm::EDGetTokenT<edm::ValueMap<float>> sigmatofpiToken_;
edm::EDGetTokenT<edm::ValueMap<float>> sigmatofkToken_;
edm::EDGetTokenT<edm::ValueMap<float>> sigmatofpToken_;
edm::EDGetTokenT<reco::VertexCollection> vtxsToken_;
edm::EDGetTokenT<edm::ValueMap<float>> trackMTDTimeQualityToken_;
const double vtxMaxSigmaT_;
Expand All @@ -71,6 +74,9 @@ TOFPIDProducer::TOFPIDProducer(const ParameterSet& iConfig)
sigmatmtdToken_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmatmtdSrc"))),
tofkToken_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofkSrc"))),
tofpToken_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofpSrc"))),
sigmatofpiToken_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmatofpiSrc"))),
sigmatofkToken_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmatofkSrc"))),
sigmatofpToken_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmatofpSrc"))),
vtxsToken_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vtxsSrc"))),
trackMTDTimeQualityToken_(
consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("trackMTDTimeQualityVMapTag"))),
Expand Down Expand Up @@ -110,6 +116,12 @@ void TOFPIDProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptio
->setComment("Input ValueMap for track tof as kaon");
desc.add<edm::InputTag>("tofpSrc", edm::InputTag("trackExtenderWithMTD:generalTrackTofP"))
->setComment("Input ValueMap for track tof as proton");
desc.add<edm::InputTag>("sigmatofpiSrc", edm::InputTag("trackExtenderWithMTD:generalTrackSigmaTofPi"))
->setComment("Input ValueMap for track sigma(tof) as pion");
desc.add<edm::InputTag>("sigmatofkSrc", edm::InputTag("trackExtenderWithMTD:generalTrackSigmaTofK"))
->setComment("Input ValueMap for track sigma(tof) as kaon");
desc.add<edm::InputTag>("sigmatofpSrc", edm::InputTag("trackExtenderWithMTD:generalTrackSigmaTofP"))
->setComment("Input ValueMap for track sigma(tof) as proton");
desc.add<edm::InputTag>("vtxsSrc", edm::InputTag("unsortedOfflinePrimaryVertices4DwithPID"))
->setComment("Input primary vertex collection");
desc.add<edm::InputTag>("trackMTDTimeQualityVMapTag", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA"))
Expand Down Expand Up @@ -163,6 +175,12 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) {

const auto& tofpIn = ev.get(tofpToken_);

const auto& sigmatofpiIn = ev.get(sigmatofpiToken_);

const auto& sigmatofkIn = ev.get(sigmatofkToken_);

const auto& sigmatofpIn = ev.get(sigmatofpToken_);

const auto& vtxs = ev.get(vtxsToken_);

const auto& trackMVAQualIn = ev.get(trackMTDTimeQualityToken_);
Expand All @@ -185,6 +203,9 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) {
float sigmat0safe = sigmat0In[trackref];
float sigmatmtd = (sigmatmtdIn[trackref] > 0. && fixedT0Error_ > 0.) ? fixedT0Error_ : sigmatmtdIn[trackref];
float sigmat0 = sigmatmtd;
float sigmatofpi = sigmatofpiIn[trackref];
float sigmatofk = sigmatofkIn[trackref];
float sigmatofp = sigmatofpIn[trackref];

float prob_pi = -1.;
float prob_k = -1.;
Expand All @@ -194,7 +215,9 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) {

if (sigmat0 > 0. && (!MVASel_ || (MVASel_ && trackMVAQual >= minTrackTimeQuality_))) {
double rsigmazsq = 1. / track.dzError() / track.dzError();
double rsigmat = 1. / sigmatmtd;
std::array<double, 3> rsigmat = {{1. / std::sqrt(sigmatmtd * sigmatmtd + sigmatofpi * sigmatofpi),
1. / std::sqrt(sigmatmtd * sigmatmtd + sigmatofk * sigmatofk),
1. / std::sqrt(sigmatmtd * sigmatmtd + sigmatofp * sigmatofp)}};

//find associated vertex
int vtxidx = -1;
Expand All @@ -217,7 +240,7 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) {
}
if (vtx.tError() > 0. && vtx.tError() < vtxMaxSigmaT_) {
double dt = std::abs(t0 - vtx.t());
double dtsig = dt * rsigmat;
double dtsig = dt * rsigmat[0]; //pion hp. uncertainty
double chisq = dz * dz * rsigmazsq + dtsig * dtsig;
if (dz < maxDz_ && dtsig < maxDtSignificance_ && chisq < minchisq) {
minchisq = chisq;
Expand Down Expand Up @@ -245,15 +268,15 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) {
const reco::Vertex& vtxnom = vtxs[vtxidx];
double dznom = std::abs(track.dz(vtxnom.position()));
double dtnom = std::abs(t0 - vtxnom.t());
double dtsignom = dtnom * rsigmat;
double dtsignom = dtnom * rsigmat[0];
double chisqnom = dznom * dznom * rsigmazsq + dtsignom * dtsignom;

//recompute t0 for alternate mass hypotheses
double t0_best = t0;

//reliable match, revert to raw mtd time uncertainty
//reliable match, revert to raw mtd time uncertainty + tof uncertainty for pion hp
if (dtsignom < maxDtSignificance_) {
sigmat0safe = sigmatmtd;
sigmat0safe = 1. / rsigmat[0];
}

double tmtd = tmtdIn[trackref];
Expand Down Expand Up @@ -284,15 +307,15 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) {
double chisqdz = dz * dz * rsigmazsq;

double dt_k = std::abs(t0_k - vtx.t());
double dtsig_k = dt_k * rsigmat;
double dtsig_k = dt_k * rsigmat[1];
double chisq_k = chisqdz + dtsig_k * dtsig_k;

if (dtsig_k < maxDtSignificance_ && chisq_k < chisqmin_k) {
chisqmin_k = chisq_k;
}

double dt_p = std::abs(t0_p - vtx.t());
double dtsig_p = dt_p * rsigmat;
double dtsig_p = dt_p * rsigmat[2];
double chisq_p = chisqdz + dtsig_p * dtsig_p;

if (dtsig_p < maxDtSignificance_ && chisq_p < chisqmin_p) {
Expand All @@ -303,13 +326,13 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) {
chisqmin = chisq_k;
t0_best = t0_k;
t0safe = t0_k;
sigmat0safe = sigmatmtd;
sigmat0safe = 1. / rsigmat[1];
}
if (dtsig_p < maxDtSignificance_ && chisq_p < chisqmin) {
chisqmin = chisq_p;
t0_best = t0_p;
t0safe = t0_p;
sigmat0safe = sigmatmtd;
sigmat0safe = 1. / rsigmat[2];
}
}

Expand Down
35 changes: 23 additions & 12 deletions RecoMTD/TrackExtender/plugins/TrackExtenderWithMTD.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -200,6 +200,7 @@ namespace {

float dt;
float dterror;
float dterror2;
float dtchi2;

float dt_best;
Expand Down Expand Up @@ -227,7 +228,7 @@ namespace {
};

enum class TofCalc { kCost = 1, kSegm = 2, kMixd = 3 };
enum class SigmaTofCalc { kCost = 1, kSegm = 2 };
enum class SigmaTofCalc { kCost = 1, kSegm = 2, kMixd = 3 };

const TrackTofPidInfo computeTrackTofPidInfo(float magp2,
float length,
Expand Down Expand Up @@ -279,6 +280,11 @@ namespace {
case SigmaTofCalc::kSegm:
res = trs.computeSigmaTof(mass_inv2);
break;
case SigmaTofCalc::kMixd:
float res1 = tofpid.pathlength * c_inv * trs.segmentSigmaMom_[trs.nSegment_ - 1] /
(magp2 * sqrt(magp2 + 1 / mass_inv2) * mass_inv2);
float res2 = trs.computeSigmaTof(mass_inv2);
res = sqrt(res1 * res1 + res2 * res2 + 2 * res1 * res2);
}

return res;
Expand All @@ -300,15 +306,18 @@ namespace {
tofpid.sigma_dt_p = sigmadeltat(m_p_inv2);

tofpid.dt = tofpid.tmtd - tofpid.dt_pi - t_vtx; //assume by default the pi hypothesis
tofpid.dterror = sqrt(tofpid.tmtderror * tofpid.tmtderror + t_vtx_err * t_vtx_err);
tofpid.dterror2 = tofpid.tmtderror * tofpid.tmtderror + t_vtx_err * t_vtx_err;
tofpid.betaerror = 0.f;
if (addPIDError) {
tofpid.dterror =
sqrt(tofpid.dterror * tofpid.dterror + (tofpid.dt_p - tofpid.dt_pi) * (tofpid.dt_p - tofpid.dt_pi));
tofpid.dterror2 = tofpid.dterror2 + (tofpid.dt_p - tofpid.dt_pi) * (tofpid.dt_p - tofpid.dt_pi);
tofpid.betaerror = tofpid.beta_p - tofpid.beta_pi;
} else {
// only add sigma(TOF) if not considering mass hp. uncertainty
tofpid.dterror2 = tofpid.dterror2 + tofpid.sigma_dt_pi * tofpid.sigma_dt_pi;
}
tofpid.dterror = sqrt(tofpid.dterror2);

tofpid.dtchi2 = (tofpid.dt * tofpid.dt) / (tofpid.dterror * tofpid.dterror);
tofpid.dtchi2 = (tofpid.dt * tofpid.dt) / tofpid.dterror2;

tofpid.dt_best = tofpid.dt;
tofpid.dterror_best = tofpid.dterror;
Expand All @@ -320,11 +329,12 @@ namespace {

if (!addPIDError) {
//*TODO* deal with heavier nucleons and/or BSM case here?
const float dterror2_wo_sigmatof = tofpid.dterror2 - tofpid.sigma_dt_pi * tofpid.sigma_dt_pi;
float chi2_pi = tofpid.dtchi2;
float chi2_k =
(tofpid.tmtd - tofpid.dt_k - t_vtx) * (tofpid.tmtd - tofpid.dt_k - t_vtx) / (tofpid.dterror * tofpid.dterror);
float chi2_p =
(tofpid.tmtd - tofpid.dt_p - t_vtx) * (tofpid.tmtd - tofpid.dt_p - t_vtx) / (tofpid.dterror * tofpid.dterror);
float chi2_k = (tofpid.tmtd - tofpid.dt_k - t_vtx) * (tofpid.tmtd - tofpid.dt_k - t_vtx) /
rappoccio marked this conversation as resolved.
Show resolved Hide resolved
(dterror2_wo_sigmatof + tofpid.sigma_dt_k * tofpid.sigma_dt_k);
float chi2_p = (tofpid.tmtd - tofpid.dt_p - t_vtx) * (tofpid.tmtd - tofpid.dt_p - t_vtx) /
(dterror2_wo_sigmatof + tofpid.sigma_dt_p * tofpid.sigma_dt_p);

float rawprob_pi = exp(-0.5f * chi2_pi);
float rawprob_k = exp(-0.5f * chi2_k);
Expand Down Expand Up @@ -1107,7 +1117,8 @@ namespace {
t_vtx,
t_vtx_err, //put vtx error by hand for the moment
false,
TofCalc::kMixd);
TofCalc::kMixd,
SigmaTofCalc::kMixd);
MTDHitMatchingInfo mi;
mi.hit = &hit;
mi.estChi2 = est.second;
Expand Down Expand Up @@ -1383,7 +1394,7 @@ reco::Track TrackExtenderWithMTDT<TrackCollection>::buildTrack(const reco::Track
//
// Protect against incompatible times
//
float err1 = tofInfo.dterror * tofInfo.dterror;
float err1 = tofInfo.dterror2;
float err2 = mtdhit2->timeError() * mtdhit2->timeError();
if (cms_rounding::roundIfNear0(err1) == 0.f || cms_rounding::roundIfNear0(err2) == 0.f) {
edm::LogError("TrackExtenderWithMTD")
Expand Down Expand Up @@ -1433,7 +1444,7 @@ reco::Track TrackExtenderWithMTDT<TrackCollection>::buildTrack(const reco::Track
tmtdOut = thit;
sigmatmtdOut = thiterror;
t0 = tofInfo.dt;
covt0t0 = tofInfo.dterror * tofInfo.dterror;
covt0t0 = tofInfo.dterror2;
betaOut = tofInfo.beta_pi;
covbetabeta = tofInfo.betaerror * tofInfo.betaerror;
tofpi = tofInfo.dt_pi;
Expand Down
8 changes: 7 additions & 1 deletion RecoVertex/PrimaryVertexProducer/interface/VertexTimeAlgorithmFromTracksPID.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ class VertexTimeAlgorithmFromTracksPID : public VertexTimeAlgorithmBase {
protected:
struct TrackInfo {
double trkWeight;
double trkTimeError;
double trkTimeErrorHyp[3];
double trkTimeHyp[3];
};

Expand All @@ -31,6 +31,9 @@ class VertexTimeAlgorithmFromTracksPID : public VertexTimeAlgorithmBase {
edm::EDGetTokenT<edm::ValueMap<float>> const trackMTDTofPiToken_;
edm::EDGetTokenT<edm::ValueMap<float>> const trackMTDTofKToken_;
edm::EDGetTokenT<edm::ValueMap<float>> const trackMTDTofPToken_;
edm::EDGetTokenT<edm::ValueMap<float>> const trackMTDSigmaTofPiToken_;
edm::EDGetTokenT<edm::ValueMap<float>> const trackMTDSigmaTofKToken_;
edm::EDGetTokenT<edm::ValueMap<float>> const trackMTDSigmaTofPToken_;

double const minTrackVtxWeight_;
double const minTrackTimeQuality_;
Expand All @@ -46,6 +49,9 @@ class VertexTimeAlgorithmFromTracksPID : public VertexTimeAlgorithmBase {
edm::ValueMap<float> trackMTDTofPi_;
edm::ValueMap<float> trackMTDTofK_;
edm::ValueMap<float> trackMTDTofP_;
edm::ValueMap<float> trackMTDSigmaTofPi_;
edm::ValueMap<float> trackMTDSigmaTofK_;
edm::ValueMap<float> trackMTDSigmaTofP_;
};

#endif
49 changes: 38 additions & 11 deletions RecoVertex/PrimaryVertexProducer/src/VertexTimeAlgorithmFromTracksPID.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,9 @@ VertexTimeAlgorithmFromTracksPID::VertexTimeAlgorithmFromTracksPID(edm::Paramete
trackMTDTofPiToken_(iCC.consumes(iConfig.getParameter<edm::InputTag>("trackMTDTofPiVMapTag"))),
trackMTDTofKToken_(iCC.consumes(iConfig.getParameter<edm::InputTag>("trackMTDTofKVMapTag"))),
trackMTDTofPToken_(iCC.consumes(iConfig.getParameter<edm::InputTag>("trackMTDTofPVMapTag"))),
trackMTDSigmaTofPiToken_(iCC.consumes(iConfig.getParameter<edm::InputTag>("trackMTDSigmaTofPiVMapTag"))),
trackMTDSigmaTofKToken_(iCC.consumes(iConfig.getParameter<edm::InputTag>("trackMTDSigmaTofKVMapTag"))),
trackMTDSigmaTofPToken_(iCC.consumes(iConfig.getParameter<edm::InputTag>("trackMTDSigmaTofPVMapTag"))),
minTrackVtxWeight_(iConfig.getParameter<double>("minTrackVtxWeight")),
minTrackTimeQuality_(iConfig.getParameter<double>("minTrackTimeQuality")),
probPion_(iConfig.getParameter<double>("probPion")),
Expand All @@ -46,6 +49,12 @@ void VertexTimeAlgorithmFromTracksPID::fillPSetDescription(edm::ParameterSetDesc
->setComment("Input ValueMap for track tof as kaon");
iDesc.add<edm::InputTag>("trackMTDTofPVMapTag", edm::InputTag("trackExtenderWithMTD:generalTrackTofP"))
->setComment("Input ValueMap for track tof as proton");
iDesc.add<edm::InputTag>("trackMTDSigmaTofPiVMapTag", edm::InputTag("trackExtenderWithMTD:generalTrackSigmaTofPi"))
->setComment("Input ValueMap for track tof uncertainty as pion");
iDesc.add<edm::InputTag>("trackMTDSigmaTofKVMapTag", edm::InputTag("trackExtenderWithMTD:generalTrackSigmaTofK"))
->setComment("Input ValueMap for track tof uncertainty as kaon");
iDesc.add<edm::InputTag>("trackMTDSigmaTofPVMapTag", edm::InputTag("trackExtenderWithMTD:generalTrackSigmaTofP"))
->setComment("Input ValueMap for track tof uncertainty as proton");

iDesc.add<double>("minTrackVtxWeight", 0.5)->setComment("Minimum track weight");
iDesc.add<double>("minTrackTimeQuality", 0.8)->setComment("Minimum MVA Quality selection on tracks");
Expand All @@ -66,6 +75,9 @@ void VertexTimeAlgorithmFromTracksPID::setEvent(edm::Event& iEvent, edm::EventSe
trackMTDTofPi_ = iEvent.get(trackMTDTofPiToken_);
trackMTDTofK_ = iEvent.get(trackMTDTofKToken_);
trackMTDTofP_ = iEvent.get(trackMTDTofPToken_);
trackMTDSigmaTofPi_ = iEvent.get(trackMTDSigmaTofPiToken_);
trackMTDSigmaTofK_ = iEvent.get(trackMTDSigmaTofKToken_);
trackMTDSigmaTofP_ = iEvent.get(trackMTDSigmaTofPToken_);
}

bool VertexTimeAlgorithmFromTracksPID::vertexTime(float& vtxTime,
Expand Down Expand Up @@ -102,23 +114,36 @@ bool VertexTimeAlgorithmFromTracksPID::vertexTime(float& vtxTime,
auto& trkInfo = v_trackInfo.back();

trkInfo.trkWeight = trkWeight;
trkInfo.trkTimeError = trkTimeError;
trkInfo.trkTimeErrorHyp[0] =
std::sqrt(trkTimeError * trkTimeError +
trackMTDSigmaTofPi_[trk.trackBaseRef()] * trackMTDSigmaTofPi_[trk.trackBaseRef()]);
trkInfo.trkTimeErrorHyp[1] =
std::sqrt(trkTimeError * trkTimeError +
trackMTDSigmaTofK_[trk.trackBaseRef()] * trackMTDSigmaTofK_[trk.trackBaseRef()]);
trkInfo.trkTimeErrorHyp[2] =
std::sqrt(trkTimeError * trkTimeError +
trackMTDSigmaTofP_[trk.trackBaseRef()] * trackMTDSigmaTofP_[trk.trackBaseRef()]);

trkInfo.trkTimeHyp[0] = trkTime - trackMTDTofPi_[trk.trackBaseRef()];
trkInfo.trkTimeHyp[1] = trkTime - trackMTDTofK_[trk.trackBaseRef()];
trkInfo.trkTimeHyp[2] = trkTime - trackMTDTofP_[trk.trackBaseRef()];

auto const wgt = trkWeight / (trkTimeError * trkTimeError);
wsum += wgt;
double const wgt[3] = {trkWeight / (trkInfo.trkTimeErrorHyp[0] * trkInfo.trkTimeErrorHyp[0]),
trkWeight / (trkInfo.trkTimeErrorHyp[1] * trkInfo.trkTimeErrorHyp[1]),
trkWeight / (trkInfo.trkTimeErrorHyp[2] * trkInfo.trkTimeErrorHyp[2])};

for (uint j = 0; j < 3; ++j) {
tsum += wgt * trkInfo.trkTimeHyp[j] * a[j];
wsum += wgt[j] * a[j];
tsum += wgt[j] * a[j] * trkInfo.trkTimeHyp[j];
}

LOG << "vertexTimeFromTracks: track"
<< " pt=" << trk.track().pt() << " eta=" << trk.track().eta() << " phi=" << trk.track().phi()
<< " vtxWeight=" << trkWeight << " time=" << trkTime << " timeError=" << trkTimeError
<< " timeQuality=" << trkTimeQuality << " timeHyp[pion]=" << trkInfo.trkTimeHyp[0]
<< " timeHyp[kaon]=" << trkInfo.trkTimeHyp[1] << " timeHyp[proton]=" << trkInfo.trkTimeHyp[2];
<< " timeQuality=" << trkTimeQuality << " timeHyp[pion]=" << trkInfo.trkTimeHyp[0] << " +/- "
<< trkInfo.trkTimeErrorHyp[0] << " timeHyp[kaon]=" << trkInfo.trkTimeHyp[1] << " +/- "
<< trkInfo.trkTimeErrorHyp[1] << " timeHyp[proton]=" << trkInfo.trkTimeHyp[2] << " +/- "
<< trkInfo.trkTimeErrorHyp[2];
}
}
}
Expand All @@ -132,27 +157,29 @@ bool VertexTimeAlgorithmFromTracksPID::vertexTime(float& vtxTime,
w2sum = 0;

for (auto const& trkInfo : v_trackInfo) {
double dt = trkInfo.trkTimeError;
double dt[3] = {trkInfo.trkTimeErrorHyp[0], trkInfo.trkTimeErrorHyp[1], trkInfo.trkTimeErrorHyp[2]};
double e[3] = {0, 0, 0};
const double cut_off = 4.5;
double Z = vdt::fast_exp(
-beta * cut_off); // outlier rejection term Z_0 = exp(-beta * cut_off) = exp(-beta * 0.5 * 3 * 3)

for (unsigned int j = 0; j < 3; j++) {
auto const tpull = (trkInfo.trkTimeHyp[j] - t0) / dt;
auto const tpull = (trkInfo.trkTimeHyp[j] - t0) / dt[j];
e[j] = vdt::fast_exp(-0.5 * beta * tpull * tpull);
Z += a[j] * e[j];
}

double wsum_trk = 0;
double wsum_trk = 0, wsum_sigma_trk = 0;
for (uint j = 0; j < 3; j++) {
double wt = a[j] * e[j] / Z;
double w = wt * trkInfo.trkWeight / (dt * dt);
double w = wt * trkInfo.trkWeight / (dt[j] * dt[j]);
wsum_trk += w;
wsum_sigma_trk += w * dt[j];
tsum += w * trkInfo.trkTimeHyp[j];
}

wsum += wsum_trk;
w2sum += wsum_trk * wsum_trk * (dt * dt) / trkInfo.trkWeight;
w2sum += wsum_sigma_trk * wsum_sigma_trk;
}

if (wsum < 1e-10) {
Expand Down
7 changes: 1 addition & 6 deletions RecoVertex/PrimaryVertexProducer/src/VertexTimeAlgorithmLegacy4D.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,6 @@ bool VertexTimeAlgorithmLegacy4D::vertexTime(float& vtxTime, float& vtxTimeError
double sumwt = 0.;
double sumwt2 = 0.;
double sumw = 0.;
double vartime = 0.;

for (const auto& trk : vtx.originalTracks()) {
const double time = trk.timeExt();
Expand All @@ -43,12 +42,8 @@ bool VertexTimeAlgorithmLegacy4D::vertexTime(float& vtxTime, float& vtxTimeError
}

if (sumw > 0) {
double sumsq = sumwt2 - sumwt * sumwt / sumw;
double chisq = num_track > 1 ? sumsq / double(num_track - 1) : sumsq / double(num_track);
vartime = chisq / sumw;

vtxTime = sumwt / sumw;
vtxTimeError = sqrt(vartime);
vtxTimeError = 1 / sqrt(sumw);
return true;
}

Expand Down