127 constexpr
double m_k = 0.493677;
128 constexpr
double m_p = 0.9382720813;
129 constexpr
double c_cm_ns = CLHEP::c_light * CLHEP::ns / CLHEP::cm;
134 const auto&
tracks = *tracksH;
138 const auto& t0In = *t0H;
142 const auto& tmtdIn = *tmtdH;
146 const auto& sigmat0In = *sigmat0H;
150 const auto& sigmatmtdIn = *sigmatmtdH;
154 const auto& pathLengthIn = *pathLengthH;
158 const auto& pIn = *pH;
162 const auto& vtxs = *vtxsH;
165 std::vector<float> t0OutRaw;
166 std::vector<float> sigmat0OutRaw;
167 std::vector<float> t0safeOutRaw;
168 std::vector<float> sigmat0safeOutRaw;
169 std::vector<float> probPiOutRaw;
170 std::vector<float> probKOutRaw;
171 std::vector<float> probPOutRaw;
174 for (
unsigned int itrack = 0; itrack <
tracks.size(); ++itrack) {
177 float t0 = t0In[trackref];
179 float sigmat0safe = sigmat0In[trackref];
181 float sigmat0 = sigmatmtd;
188 double rsigmazsq = 1. /
track.dzError() /
track.dzError();
189 double rsigmat = 1. / sigmatmtd;
193 int vtxidxmindz = -1;
194 int vtxidxminchisq = -1;
198 for (
unsigned int ivtx = 0; ivtx < vtxs.size(); ++ivtx) {
200 float w =
vtx.trackWeight(trackref);
212 double dtsig =
dt * rsigmat;
213 double chisq =
dz *
dz * rsigmazsq + dtsig * dtsig;
216 vtxidxminchisq = ivtx;
225 if (vtxidxmindz >= 0 && !(vtxs[vtxidxmindz].tError() > 0. && vtxs[vtxidxmindz].tError() <
vtxMaxSigmaT_)) {
226 vtxidx = vtxidxmindz;
227 }
else if (vtxidxminchisq >= 0) {
228 vtxidx = vtxidxminchisq;
229 }
else if (vtxidxmindz >= 0) {
230 vtxidx = vtxidxmindz;
235 if (vtxidx >= 0 && vtxs[vtxidx].tError() > 0. && vtxs[vtxidx].tError() <
vtxMaxSigmaT_) {
240 double dtsignom = dtnom * rsigmat;
241 double chisqnom = dznom * dznom * rsigmazsq + dtsignom * dtsignom;
248 sigmat0safe = sigmatmtd;
251 double tmtd = tmtdIn[trackref];
252 double pathlength = pathLengthIn[trackref];
253 double magp = pIn[trackref];
255 double gammasq_k = 1. + magp * magp / m_k / m_k;
256 double beta_k =
std::sqrt(1. - 1. / gammasq_k);
257 double t0_k = tmtd - pathlength / beta_k *
c_inv;
259 double gammasq_p = 1. + magp * magp / m_p / m_p;
260 double beta_p =
std::sqrt(1. - 1. / gammasq_p);
261 double t0_p = tmtd - pathlength / beta_p *
c_inv;
263 double chisqmin = chisqnom;
265 double chisqmin_pi = chisqnom;
279 double chisqdz =
dz *
dz * rsigmazsq;
282 double dtsig_k = dt_k * rsigmat;
283 double chisq_k = chisqdz + dtsig_k * dtsig_k;
286 chisqmin_k = chisq_k;
290 double dtsig_p = dt_p * rsigmat;
291 double chisq_p = chisqdz + dtsig_p * dtsig_p;
294 chisqmin_p = chisq_p;
301 sigmat0safe = sigmatmtd;
307 sigmat0safe = sigmatmtd;
313 double rawprob_pi =
exp(-0.5 * chisqmin_pi);
314 double rawprob_k =
exp(-0.5 * chisqmin_k);
315 double rawprob_p =
exp(-0.5 * chisqmin_p);
317 double normprob = 1. / (rawprob_pi + rawprob_k + rawprob_p);
319 prob_pi = rawprob_pi * normprob;
320 prob_k = rawprob_k * normprob;
321 prob_p = rawprob_p * normprob;
323 double prob_heavy = 1. - prob_pi;
331 t0OutRaw.push_back(
t0);
332 sigmat0OutRaw.push_back(sigmat0);
333 t0safeOutRaw.push_back(t0safe);
334 sigmat0safeOutRaw.push_back(sigmat0safe);
335 probPiOutRaw.push_back(prob_pi);
336 probKOutRaw.push_back(prob_k);
337 probPOutRaw.push_back(prob_p);