Skip to content

Commit fdd6e31

Browse files
noepalmnoepalm
committed
Added sigma(p) to TrackSegments; added sigma(TOF) computation, now saved to event
1 parent 7d488e0 commit fdd6e31

File tree

1 file changed

+139
-18
lines changed

1 file changed

+139
-18
lines changed

‎RecoMTD/TrackExtender/plugins/TrackExtenderWithMTD.cc

Lines changed: 139 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -98,13 +98,16 @@ namespace {
9898
public:
9999
TrackSegments() = default;
100100

101-
inline uint32_t addSegment(float tPath, float tMom2) {
101+
inline uint32_t addSegment(float tPath, float tMom2, float sigmaMom) {
102102
segmentPathOvc_.emplace_back(tPath * c_inv);
103103
segmentMom2_.emplace_back(tMom2);
104+
segmentSigmaMom_.emplace_back(sigmaMom);
104105
nSegment_++;
105106

106107
LogTrace("TrackExtenderWithMTD") << "addSegment # " << nSegment_ << " s = " << tPath
107-
<< " p = " << std::sqrt(tMom2);
108+
<< " p = " << std::sqrt(tMom2)
109+
<< " sigma_p = " << sigmaMom
110+
<< " sigma_p/p = " << sigmaMom / std::sqrt(tMom2) * 100 << " %";
108111

109112
return nSegment_;
110113
}
@@ -116,12 +119,28 @@ namespace {
116119
float beta = std::sqrt(1.f - 1.f / gammasq);
117120
tof += segmentPathOvc_[iSeg] / beta;
118121

119-
LogTrace("TrackExtenderWithMTD") << " TOF Segment # " << iSeg + 1 << " p = " << std::sqrt(segmentMom2_[iSeg])
120-
<< " tof = " << tof;
122+
float sigma_tof = segmentPathOvc_[iSeg] * segmentSigmaMom_[iSeg] / (segmentMom2_[iSeg] * sqrt(segmentMom2_[iSeg] + 1/mass_inv2) * mass_inv2);
123+
LogTrace("TrackExtenderWithMTD") << "TOF Segment # " << iSeg + 1 << " p = " << std::sqrt(segmentMom2_[iSeg]) << "\t1/gamma_sq = " << 1/gammasq
124+
<< "\tdelta(tof) = " << segmentPathOvc_[iSeg] / beta << "\tsigma_delta(tof) = " << sigma_tof << "\tsigma_tof/delta(tof) = " << sigma_tof / (segmentPathOvc_[iSeg] / beta) * 100 << "%\ttof = " << tof;
121125
}
122126

123127
return tof;
124128
}
129+
130+
inline float computeSigmaTof(float mass_inv2) const {
131+
float sigmatof = 0.;
132+
for (uint32_t iSeg = 0; iSeg < nSegment_; iSeg++) {
133+
134+
float sigma_i = segmentPathOvc_[iSeg] * segmentSigmaMom_[iSeg] / (segmentMom2_[iSeg] * sqrt(segmentMom2_[iSeg] + 1/mass_inv2) * mass_inv2);
135+
136+
for (uint32_t jSeg = 0; jSeg < nSegment_; jSeg++) {
137+
float sigma_j = segmentPathOvc_[jSeg] * segmentSigmaMom_[jSeg] / (segmentMom2_[jSeg] * sqrt(segmentMom2_[jSeg] + 1/mass_inv2) * mass_inv2);
138+
sigmatof += sigma_i * sigma_j;
139+
}
140+
}
141+
142+
return sqrt(sigmatof);
143+
}
125144

126145
inline uint32_t size() const { return nSegment_; }
127146

@@ -144,6 +163,7 @@ namespace {
144163
uint32_t nSegment_ = 0;
145164
std::vector<float> segmentPathOvc_;
146165
std::vector<float> segmentMom2_;
166+
std::vector<float> segmentSigmaMom_;
147167
};
148168

149169
struct TrackTofPidInfo {
@@ -164,21 +184,25 @@ namespace {
164184
float gammasq_pi;
165185
float beta_pi;
166186
float dt_pi;
187+
float sigma_dt_pi;
167188

168189
float gammasq_k;
169190
float beta_k;
170191
float dt_k;
192+
float sigma_dt_k;
171193

172194
float gammasq_p;
173195
float beta_p;
174196
float dt_p;
197+
float sigma_dt_p;
175198

176199
float prob_pi;
177200
float prob_k;
178201
float prob_p;
179202
};
180203

181204
enum class TofCalc { kCost = 1, kSegm = 2, kMixd = 3 };
205+
enum class SigmaTofCalc { kCost = 1, kSegm = 2 };
182206

183207
const TrackTofPidInfo computeTrackTofPidInfo(float magp2,
184208
float length,
@@ -188,7 +212,8 @@ namespace {
188212
float t_vtx,
189213
float t_vtx_err,
190214
bool addPIDError = true,
191-
TofCalc choice = TofCalc::kCost) {
215+
TofCalc choice = TofCalc::kCost,
216+
SigmaTofCalc sigma_choice = SigmaTofCalc::kCost) {
192217
constexpr float m_pi = 0.13957018f;
193218
constexpr float m_pi_inv2 = 1.0f / m_pi / m_pi;
194219
constexpr float m_k = 0.493677f;
@@ -218,17 +243,38 @@ namespace {
218243
return res;
219244
};
220245

246+
auto sigmadeltat = [&](const float mass_inv2) {
247+
float res(1.f);
248+
switch (sigma_choice) {
249+
case SigmaTofCalc::kCost:
250+
// sigma(t) = sigma(p) * |dt/dp| = sigma(p) * DeltaL/c * m^2 / (p^2 * E)
251+
res = tofpid.pathlength * c_inv * trs.segmentSigmaMom_[trs.nSegment_ - 1] / (magp2 * sqrt(magp2 + 1/mass_inv2) * mass_inv2);
252+
break;
253+
case SigmaTofCalc::kSegm:
254+
res = trs.computeSigmaTof(mass_inv2);
255+
break;
256+
}
257+
258+
return res;
259+
};
260+
261+
LogTrace("TrackExtenderWithMTD") << "Pion:";
221262
tofpid.gammasq_pi = 1.f + magp2 * m_pi_inv2;
222263
tofpid.beta_pi = std::sqrt(1.f - 1.f / tofpid.gammasq_pi);
223264
tofpid.dt_pi = deltat(m_pi_inv2, tofpid.beta_pi);
265+
tofpid.sigma_dt_pi = sigmadeltat(m_pi_inv2);
224266

267+
LogTrace("TrackExtenderWithMTD") << "Kaon:";
225268
tofpid.gammasq_k = 1.f + magp2 * m_k_inv2;
226269
tofpid.beta_k = std::sqrt(1.f - 1.f / tofpid.gammasq_k);
227270
tofpid.dt_k = deltat(m_k_inv2, tofpid.beta_k);
271+
tofpid.sigma_dt_k = sigmadeltat(m_k_inv2);
228272

273+
LogTrace("TrackExtenderWithMTD") << "Proton:";
229274
tofpid.gammasq_p = 1.f + magp2 * m_p_inv2;
230275
tofpid.beta_p = std::sqrt(1.f - 1.f / tofpid.gammasq_p);
231276
tofpid.dt_p = deltat(m_p_inv2, tofpid.beta_p);
277+
tofpid.sigma_dt_p = sigmadeltat(m_p_inv2);
232278

233279
tofpid.dt = tofpid.tmtd - tofpid.dt_pi - t_vtx; //assume by default the pi hypothesis
234280
tofpid.dterror = sqrt(tofpid.tmtderror * tofpid.tmtderror + t_vtx_err * t_vtx_err);
@@ -323,7 +369,28 @@ namespace {
323369
validpropagation = false;
324370
}
325371
pathlength1 += layerpathlength;
326-
trs.addSegment(layerpathlength, (it + 1)->updatedState().globalMomentum().mag2());
372+
373+
// // sigma(p) using cartesian error
374+
// float sigma_p_cartesian = 0;
375+
// float p[3] = {(it + 1)->updatedState().globalMomentum().x(), (it + 1)->updatedState().globalMomentum().y(),
376+
// (it + 1)->updatedState().globalMomentum().z()};
377+
// // calculate sigma_p on p = sqrt(p_x^2 + p_y^2 + p_z^2) as:
378+
// // sigma_p = sqrt(sigma_px^2 * p_x^2 + [same for y, z] + cov(px, py) * px * py + [same for x-z, y-z]) / p
379+
// for(int i = 3; i < 6; i++){
380+
// for(int j = 3; j < 6; j++){
381+
// // summing 2 * sigma_pij * p_i * p_j / p^2
382+
// // std::cout << "sigma_p_" << i << j << " = " << (it + 1)->updatedState().cartesianError().matrix()(i, j) << " " << "p_" << i << " = " << p[i-3] << " " << "p_" << j << " = " << p[j-3] << " ";
383+
// sigma_p_cartesian += (it + 1)->updatedState().cartesianError().matrix()(i, j) * p[i-3] * p[j-3] / (it + 1)->updatedState().globalMomentum().mag2();
384+
// }
385+
// }
386+
387+
// sigma_p_cartesian = sqrt(sigma_p_cartesian);
388+
389+
// sigma(p) using curvilinear error (on q/p)
390+
float sigma_p = sqrt((it + 1)->updatedState().curvilinearError().matrix()(0, 0)) * (it + 1)->updatedState().globalMomentum().mag2();
391+
392+
trs.addSegment(layerpathlength, (it + 1)->updatedState().globalMomentum().mag2(), sigma_p);
393+
327394
LogTrace("TrackExtenderWithMTD") << "TSOS " << std::fixed << std::setw(4) << trs.size() << " R_i " << std::fixed
328395
<< std::setw(14) << it->updatedState().globalPosition().perp() << " z_i "
329396
<< std::fixed << std::setw(14) << it->updatedState().globalPosition().z()
@@ -332,7 +399,14 @@ namespace {
332399
<< std::setw(14) << (it + 1)->updatedState().globalPosition().z() << " p "
333400
<< std::fixed << std::setw(14) << (it + 1)->updatedState().globalMomentum().mag()
334401
<< " dp " << std::fixed << std::setw(14)
335-
<< (it + 1)->updatedState().globalMomentum().mag() - oldp;
402+
<< (it + 1)->updatedState().globalMomentum().mag() - oldp
403+
<< " sigma_p (from curvilinear) = " << std::fixed << std::setw(14)
404+
<< sigma_p
405+
// << " sigma_p (cartesian) = " << std::fixed << std::setw(14)
406+
// << sigma_p_cartesian
407+
<< " sigma_p/p = " << std::fixed << std::setw(14)
408+
<< sigma_p / (it + 1)->updatedState().globalMomentum().mag() * 100 << " %";
409+
336410
oldp = (it + 1)->updatedState().globalMomentum().mag();
337411
}
338412

@@ -345,12 +419,21 @@ namespace {
345419
validpropagation = false;
346420
}
347421
pathlength = pathlength1 + pathlength2;
348-
trs.addSegment(pathlength2, tscblPCA.momentum().mag2());
422+
423+
float sigma_p = sqrt(tscblPCA.curvilinearError().matrix()(0, 0)) * tscblPCA.momentum().mag2();
424+
425+
trs.addSegment(pathlength2, tscblPCA.momentum().mag2(), sigma_p);
426+
349427
LogTrace("TrackExtenderWithMTD") << "TSOS " << std::fixed << std::setw(4) << trs.size() << " R_e " << std::fixed
350428
<< std::setw(14) << tscblPCA.position().perp() << " z_e " << std::fixed
351429
<< std::setw(14) << tscblPCA.position().z() << " p " << std::fixed << std::setw(14)
352430
<< tscblPCA.momentum().mag() << " dp " << std::fixed << std::setw(14)
353-
<< tscblPCA.momentum().mag() - oldp;
431+
<< tscblPCA.momentum().mag() - oldp
432+
<< " sigma_p = " << std::fixed << std::setw(14)
433+
<< sigma_p
434+
<< " sigma_p/p = " << std::fixed << std::setw(14)
435+
<< sigma_p / tscblPCA.momentum().mag() * 100 << " %";
436+
354437
return validpropagation;
355438
}
356439

@@ -459,7 +542,10 @@ class TrackExtenderWithMTDT : public edm::stream::EDProducer<> {
459542
float& sigmatmtdOut,
460543
float& tofpi,
461544
float& tofk,
462-
float& tofp) const;
545+
float& tofp,
546+
float& sigmatofpi,
547+
float& sigmatofk,
548+
float& sigmatofp) const;
463549
reco::TrackExtra buildTrackExtra(const Trajectory& trajectory) const;
464550

465551
string dumpLayer(const DetLayer* layer) const;
@@ -481,6 +567,9 @@ class TrackExtenderWithMTDT : public edm::stream::EDProducer<> {
481567
edm::EDPutToken tofpiOrigTrkToken_;
482568
edm::EDPutToken tofkOrigTrkToken_;
483569
edm::EDPutToken tofpOrigTrkToken_;
570+
edm::EDPutToken sigmatofpiOrigTrkToken_;
571+
edm::EDPutToken sigmatofkOrigTrkToken_;
572+
edm::EDPutToken sigmatofpOrigTrkToken_;
484573
edm::EDPutToken assocOrigTrkToken_;
485574

486575
edm::EDGetTokenT<InputCollection> tracksToken_;
@@ -569,6 +658,9 @@ TrackExtenderWithMTDT<TrackCollection>::TrackExtenderWithMTDT(const ParameterSet
569658
tofpiOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackTofPi");
570659
tofkOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackTofK");
571660
tofpOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackTofP");
661+
sigmatofpiOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackSigmaTofPi");
662+
sigmatofkOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackSigmaTofK");
663+
sigmatofpOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackSigmaTofP");
572664
assocOrigTrkToken_ = produces<edm::ValueMap<int>>("generalTrackassoc");
573665

574666
builderToken_ = esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", transientTrackBuilder_));
@@ -683,6 +775,9 @@ void TrackExtenderWithMTDT<TrackCollection>::produce(edm::Event& ev, const edm::
683775
std::vector<float> tofpiOrigTrkRaw;
684776
std::vector<float> tofkOrigTrkRaw;
685777
std::vector<float> tofpOrigTrkRaw;
778+
std::vector<float> sigmatofpiOrigTrkRaw;
779+
std::vector<float> sigmatofkOrigTrkRaw;
780+
std::vector<float> sigmatofpOrigTrkRaw;
686781
std::vector<int> assocOrigTrkRaw;
687782

688783
auto const tracksH = ev.getHandle(tracksToken_);
@@ -727,6 +822,9 @@ void TrackExtenderWithMTDT<TrackCollection>::produce(edm::Event& ev, const edm::
727822

728823
LogTrace("TrackExtenderWithMTD") << "TrackExtenderWithMTD: extrapolating track " << itrack
729824
<< " p/pT = " << track->p() << " " << track->pt() << " eta = " << track->eta();
825+
LogTrace("TrackExtenderWithMTD") << "TrackExtenderWithMTD: sigma_p = "
826+
<< sqrt(track->covariance()(0,0)) * track->p2()
827+
<< " sigma_p/p = " << sqrt(track->covariance()(0,0)) * track->p() * 100 << " %";
730828

731829
float trackVtxTime = 0.f;
732830
if (useVertex_) {
@@ -803,12 +901,14 @@ void TrackExtenderWithMTDT<TrackCollection>::produce(edm::Event& ev, const edm::
803901
const auto& trajwithmtd =
804902
mtdthits.empty() ? std::vector<Trajectory>(1, trajs) : theTransformer->transform(ttrack, thits);
805903
float pMap = 0.f, betaMap = 0.f, t0Map = 0.f, sigmat0Map = -1.f, pathLengthMap = -1.f, tmtdMap = 0.f,
806-
sigmatmtdMap = -1.f, tofpiMap = 0.f, tofkMap = 0.f, tofpMap = 0.f;
904+
sigmatmtdMap = -1.f, tofpiMap = 0.f, tofkMap = 0.f, tofpMap = 0.f, sigmatofpiMap = -1.f,
905+
sigmatofkMap = -1.f, sigmatofpMap = -1.f;
807906
int iMap = -1;
808907

809908
for (const auto& trj : trajwithmtd) {
810909
const auto& thetrj = (updateTraj_ ? trj : trajs);
811-
float pathLength = 0.f, tmtd = 0.f, sigmatmtd = -1.f, tofpi = 0.f, tofk = 0.f, tofp = 0.f;
910+
float pathLength = 0.f, tmtd = 0.f, sigmatmtd = -1.f, tofpi = 0.f, tofk = 0.f, tofp = 0.f, sigmatofpi = -1.f,
911+
sigmatofk = -1.f, sigmatofp = -1.f;
812912
LogTrace("TrackExtenderWithMTD") << "TrackExtenderWithMTD: refit track " << itrack << " p/pT = " << track->p()
813913
<< " " << track->pt() << " eta = " << track->eta();
814914
reco::Track result = buildTrack(track,
@@ -823,7 +923,10 @@ void TrackExtenderWithMTDT<TrackCollection>::produce(edm::Event& ev, const edm::
823923
sigmatmtd,
824924
tofpi,
825925
tofk,
826-
tofp);
926+
tofp,
927+
sigmatofpi,
928+
sigmatofk,
929+
sigmatofp);
827930
if (result.ndof() >= 0) {
828931
/// setup the track extras
829932
reco::TrackExtra::TrajParams trajParams;
@@ -856,6 +959,9 @@ void TrackExtenderWithMTDT<TrackCollection>::produce(edm::Event& ev, const edm::
856959
tofpiMap = tofpi;
857960
tofkMap = tofk;
858961
tofpMap = tofp;
962+
sigmatofpiMap = sigmatofpi;
963+
sigmatofkMap = sigmatofk;
964+
sigmatofpMap = sigmatofp;
859965
reco::TrackExtraRef extraRef(extrasRefProd, extras->size() - 1);
860966
backtrack.setExtra((updateExtra_ ? extraRef : track->extra()));
861967
for (unsigned ihit = hitsstart; ihit < hitsend; ++ihit) {
@@ -864,8 +970,10 @@ void TrackExtenderWithMTDT<TrackCollection>::produce(edm::Event& ev, const edm::
864970
npixBarrel.push_back(backtrack.hitPattern().numberOfValidPixelBarrelHits());
865971
npixEndcap.push_back(backtrack.hitPattern().numberOfValidPixelEndcapHits());
866972
LogTrace("TrackExtenderWithMTD") << "TrackExtenderWithMTD: tmtd " << tmtdMap << " +/- " << sigmatmtdMap
867-
<< " t0 " << t0Map << " +/- " << sigmat0Map << " tof pi/K/p " << tofpiMap
868-
<< " " << tofkMap << " " << tofpMap;
973+
<< " t0 " << t0Map << " +/- " << sigmat0Map << " tof pi/K/p "
974+
<< tofpiMap << "+/-" << fmt::format("{:0.2g}", sigmatofpiMap) << " (" << fmt::format("{:0.2g}", sigmatofpiMap/tofpiMap*100) << "%) "
975+
<< tofkMap << "+/-" << fmt::format("{:0.2g}", sigmatofkMap) << " (" << fmt::format("{:0.2g}", sigmatofkMap/tofkMap*100) << "%) "
976+
<< tofpMap << "+/-" << fmt::format("{:0.2g}", sigmatofpMap) << " (" << fmt::format("{:0.2g}", sigmatofpMap/tofpMap*100) << "%) ";
869977
} else {
870978
LogTrace("TrackExtenderWithMTD") << "Error in the MTD track refitting. This should not happen";
871979
}
@@ -881,6 +989,9 @@ void TrackExtenderWithMTDT<TrackCollection>::produce(edm::Event& ev, const edm::
881989
tofpiOrigTrkRaw.push_back(tofpiMap);
882990
tofkOrigTrkRaw.push_back(tofkMap);
883991
tofpOrigTrkRaw.push_back(tofpMap);
992+
sigmatofpiOrigTrkRaw.push_back(sigmatofpiMap);
993+
sigmatofkOrigTrkRaw.push_back(sigmatofkMap);
994+
sigmatofpOrigTrkRaw.push_back(sigmatofpMap);
884995
assocOrigTrkRaw.push_back(iMap);
885996

886997
if (iMap == -1) {
@@ -915,6 +1026,9 @@ void TrackExtenderWithMTDT<TrackCollection>::produce(edm::Event& ev, const edm::
9151026
fillValueMap(ev, tracksH, tofpiOrigTrkRaw, tofpiOrigTrkToken_);
9161027
fillValueMap(ev, tracksH, tofkOrigTrkRaw, tofkOrigTrkToken_);
9171028
fillValueMap(ev, tracksH, tofpOrigTrkRaw, tofpOrigTrkToken_);
1029+
fillValueMap(ev, tracksH, sigmatofpiOrigTrkRaw, sigmatofpiOrigTrkToken_);
1030+
fillValueMap(ev, tracksH, sigmatofkOrigTrkRaw, sigmatofkOrigTrkToken_);
1031+
fillValueMap(ev, tracksH, sigmatofpOrigTrkRaw, sigmatofpOrigTrkToken_);
9181032
fillValueMap(ev, tracksH, assocOrigTrkRaw, assocOrigTrkToken_);
9191033
}
9201034

@@ -1176,7 +1290,11 @@ reco::Track TrackExtenderWithMTDT<TrackCollection>::buildTrack(const reco::Track
11761290
float& sigmatmtdOut,
11771291
float& tofpi,
11781292
float& tofk,
1179-
float& tofp) const {
1293+
float& tofp,
1294+
float& sigmatofpi,
1295+
float& sigmatofk,
1296+
float& sigmatofp
1297+
) const {
11801298
TrajectoryStateClosestToBeamLine tscbl;
11811299
bool tsbcl_status = getTrajectoryStateClosestToBeamLine(traj, bs, thePropagator, tscbl);
11821300

@@ -1308,7 +1426,7 @@ reco::Track TrackExtenderWithMTDT<TrackCollection>::buildTrack(const reco::Track
13081426
if (validmtd && validpropagation) {
13091427
//here add the PID uncertainty for later use in the 1st step of 4D vtx reconstruction
13101428
TrackTofPidInfo tofInfo =
1311-
computeTrackTofPidInfo(p.mag2(), pathlength, trs, thit, thiterror, 0.f, 0.f, true, TofCalc::kSegm);
1429+
computeTrackTofPidInfo(p.mag2(), pathlength, trs, thit, thiterror, 0.f, 0.f, true, TofCalc::kSegm, SigmaTofCalc::kCost);
13121430
pathLengthOut = pathlength; // set path length if we've got a timing hit
13131431
tmtdOut = thit;
13141432
sigmatmtdOut = thiterror;
@@ -1319,6 +1437,9 @@ reco::Track TrackExtenderWithMTDT<TrackCollection>::buildTrack(const reco::Track
13191437
tofpi = tofInfo.dt_pi;
13201438
tofk = tofInfo.dt_k;
13211439
tofp = tofInfo.dt_p;
1440+
sigmatofpi = tofInfo.sigma_dt_pi;
1441+
sigmatofk = tofInfo.sigma_dt_k;
1442+
sigmatofp = tofInfo.sigma_dt_p;
13221443
}
13231444
}
13241445

@@ -1426,4 +1547,4 @@ string TrackExtenderWithMTDT<TrackCollection>::dumpLayer(const DetLayer* layer)
14261547
#include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
14271548
typedef TrackExtenderWithMTDT<reco::TrackCollection> TrackExtenderWithMTD;
14281549

1429-
DEFINE_FWK_MODULE(TrackExtenderWithMTD);
1550+
DEFINE_FWK_MODULE(TrackExtenderWithMTD);

0 commit comments

Comments
 (0)