15 #include "CLHEP/Units/GlobalPhysicalConstants.h"
16 #include "CLHEP/Units/GlobalSystemOfUnits.h"
27 template <
class H,
class T>
30 const std::vector<T>& vec,
36 static constexpr
char t0Name[] =
"t0";
37 static constexpr
char sigmat0Name[] =
"sigmat0";
38 static constexpr
char t0safeName[] =
"t0safe";
39 static constexpr
char sigmat0safeName[] =
"sigmat0safe";
40 static constexpr
char probPiName[] =
"probPi";
41 static constexpr
char probKName[] =
"probK";
42 static constexpr
char probPName[] =
"probP";
68 vtxMaxSigmaT_(iConfig.getParameter<double>(
"vtxMaxSigmaT")),
69 maxDz_(iConfig.getParameter<double>(
"maxDz")),
70 maxDtSignificance_(iConfig.getParameter<double>(
"maxDtSignificance")),
71 minProbHeavy_(iConfig.getParameter<double>(
"minProbHeavy")),
72 fixedT0Error_(iConfig.getParameter<double>(
"fixedT0Error")) {
73 produces<edm::ValueMap<float>>(
t0Name);
78 produces<edm::ValueMap<float>>(
probKName);
79 produces<edm::ValueMap<float>>(
probPName);
87 ->setComment(
"Input ValueMap for track time at beamline");
89 ->setComment(
"Input ValueMap for track time at MTD");
91 ->setComment(
"Input ValueMap for track time uncertainty at beamline");
93 ->setComment(
"Input ValueMap for track time uncertainty at MTD");
95 ->setComment(
"Input ValueMap for track path lengh from beamline to MTD");
97 ->setComment(
"Input ValueMap for track momentum magnitude (normally from refit with MTD hits)");
99 ->setComment(
"Input primary vertex collection");
100 desc.add<
double>(
"vtxMaxSigmaT", 0.025)
101 ->setComment(
"Maximum primary vertex time uncertainty for use in particle id [ns]");
102 desc.add<
double>(
"maxDz", 0.1)
103 ->setComment(
"Maximum distance in z for track-primary vertex association for particle id [cm]");
104 desc.add<
double>(
"maxDtSignificance", 5.0)
106 "Maximum distance in time (normalized by uncertainty) for track-primary vertex association for particle id");
107 desc.add<
double>(
"minProbHeavy", 0.75)
108 ->setComment(
"Minimum probability for a particle to be a kaon or proton before reassigning the timestamp");
109 desc.add<
double>(
"fixedT0Error", 0.)->setComment(
"Use a fixed T0 uncertainty [ns]");
111 descriptions.
add(
"tofPIDProducer",
desc);
114 template <
class H,
class T>
117 const std::vector<T>& vec,
119 auto out = std::make_unique<edm::ValueMap<T>>();
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);