From fdd6e311a0f33ae5d49a9b025976a0d6f64d9e18 Mon Sep 17 00:00:00 2001 From: noepalm Date: Tue, 9 Jan 2024 16:08:46 +0100 Subject: [PATCH] Added sigma(p) to TrackSegments; added sigma(TOF) computation, now saved to event --- .../plugins/TrackExtenderWithMTD.cc | 157 ++++++++++++++++-- 1 file changed, 139 insertions(+), 18 deletions(-) diff --git a/RecoMTD/TrackExtender/plugins/TrackExtenderWithMTD.cc b/RecoMTD/TrackExtender/plugins/TrackExtenderWithMTD.cc index 3adf7eab66f03..2db3423b54aec 100644 --- a/RecoMTD/TrackExtender/plugins/TrackExtenderWithMTD.cc +++ b/RecoMTD/TrackExtender/plugins/TrackExtenderWithMTD.cc @@ -98,13 +98,16 @@ namespace { public: TrackSegments() = default; - inline uint32_t addSegment(float tPath, float tMom2) { + inline uint32_t addSegment(float tPath, float tMom2, float sigmaMom) { segmentPathOvc_.emplace_back(tPath * c_inv); segmentMom2_.emplace_back(tMom2); + segmentSigmaMom_.emplace_back(sigmaMom); nSegment_++; LogTrace("TrackExtenderWithMTD") << "addSegment # " << nSegment_ << " s = " << tPath - << " p = " << std::sqrt(tMom2); + << " p = " << std::sqrt(tMom2) + << " sigma_p = " << sigmaMom + << " sigma_p/p = " << sigmaMom / std::sqrt(tMom2) * 100 << " %"; return nSegment_; } @@ -116,12 +119,28 @@ namespace { float beta = std::sqrt(1.f - 1.f / gammasq); tof += segmentPathOvc_[iSeg] / beta; - LogTrace("TrackExtenderWithMTD") << " TOF Segment # " << iSeg + 1 << " p = " << std::sqrt(segmentMom2_[iSeg]) - << " tof = " << tof; + float sigma_tof = segmentPathOvc_[iSeg] * segmentSigmaMom_[iSeg] / (segmentMom2_[iSeg] * sqrt(segmentMom2_[iSeg] + 1/mass_inv2) * mass_inv2); + LogTrace("TrackExtenderWithMTD") << "TOF Segment # " << iSeg + 1 << " p = " << std::sqrt(segmentMom2_[iSeg]) << "\t1/gamma_sq = " << 1/gammasq + << "\tdelta(tof) = " << segmentPathOvc_[iSeg] / beta << "\tsigma_delta(tof) = " << sigma_tof << "\tsigma_tof/delta(tof) = " << sigma_tof / (segmentPathOvc_[iSeg] / beta) * 100 << "%\ttof = " << tof; } return tof; } + + inline float computeSigmaTof(float mass_inv2) const { + float sigmatof = 0.; + for (uint32_t iSeg = 0; iSeg < nSegment_; iSeg++) { + + float sigma_i = segmentPathOvc_[iSeg] * segmentSigmaMom_[iSeg] / (segmentMom2_[iSeg] * sqrt(segmentMom2_[iSeg] + 1/mass_inv2) * mass_inv2); + + for (uint32_t jSeg = 0; jSeg < nSegment_; jSeg++) { + float sigma_j = segmentPathOvc_[jSeg] * segmentSigmaMom_[jSeg] / (segmentMom2_[jSeg] * sqrt(segmentMom2_[jSeg] + 1/mass_inv2) * mass_inv2); + sigmatof += sigma_i * sigma_j; + } + } + + return sqrt(sigmatof); + } inline uint32_t size() const { return nSegment_; } @@ -144,6 +163,7 @@ namespace { uint32_t nSegment_ = 0; std::vector segmentPathOvc_; std::vector segmentMom2_; + std::vector segmentSigmaMom_; }; struct TrackTofPidInfo { @@ -164,14 +184,17 @@ namespace { float gammasq_pi; float beta_pi; float dt_pi; + float sigma_dt_pi; float gammasq_k; float beta_k; float dt_k; + float sigma_dt_k; float gammasq_p; float beta_p; float dt_p; + float sigma_dt_p; float prob_pi; float prob_k; @@ -179,6 +202,7 @@ namespace { }; enum class TofCalc { kCost = 1, kSegm = 2, kMixd = 3 }; + enum class SigmaTofCalc { kCost = 1, kSegm = 2 }; const TrackTofPidInfo computeTrackTofPidInfo(float magp2, float length, @@ -188,7 +212,8 @@ namespace { float t_vtx, float t_vtx_err, bool addPIDError = true, - TofCalc choice = TofCalc::kCost) { + TofCalc choice = TofCalc::kCost, + SigmaTofCalc sigma_choice = SigmaTofCalc::kCost) { constexpr float m_pi = 0.13957018f; constexpr float m_pi_inv2 = 1.0f / m_pi / m_pi; constexpr float m_k = 0.493677f; @@ -218,17 +243,38 @@ namespace { return res; }; + auto sigmadeltat = [&](const float mass_inv2) { + float res(1.f); + switch (sigma_choice) { + case SigmaTofCalc::kCost: + // sigma(t) = sigma(p) * |dt/dp| = sigma(p) * DeltaL/c * m^2 / (p^2 * E) + res = tofpid.pathlength * c_inv * trs.segmentSigmaMom_[trs.nSegment_ - 1] / (magp2 * sqrt(magp2 + 1/mass_inv2) * mass_inv2); + break; + case SigmaTofCalc::kSegm: + res = trs.computeSigmaTof(mass_inv2); + break; + } + + return res; + }; + + LogTrace("TrackExtenderWithMTD") << "Pion:"; tofpid.gammasq_pi = 1.f + magp2 * m_pi_inv2; tofpid.beta_pi = std::sqrt(1.f - 1.f / tofpid.gammasq_pi); tofpid.dt_pi = deltat(m_pi_inv2, tofpid.beta_pi); + tofpid.sigma_dt_pi = sigmadeltat(m_pi_inv2); + LogTrace("TrackExtenderWithMTD") << "Kaon:"; tofpid.gammasq_k = 1.f + magp2 * m_k_inv2; tofpid.beta_k = std::sqrt(1.f - 1.f / tofpid.gammasq_k); tofpid.dt_k = deltat(m_k_inv2, tofpid.beta_k); + tofpid.sigma_dt_k = sigmadeltat(m_k_inv2); + LogTrace("TrackExtenderWithMTD") << "Proton:"; tofpid.gammasq_p = 1.f + magp2 * m_p_inv2; tofpid.beta_p = std::sqrt(1.f - 1.f / tofpid.gammasq_p); tofpid.dt_p = deltat(m_p_inv2, tofpid.beta_p); + 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); @@ -323,7 +369,28 @@ namespace { validpropagation = false; } pathlength1 += layerpathlength; - trs.addSegment(layerpathlength, (it + 1)->updatedState().globalMomentum().mag2()); + + // // sigma(p) using cartesian error + // float sigma_p_cartesian = 0; + // float p[3] = {(it + 1)->updatedState().globalMomentum().x(), (it + 1)->updatedState().globalMomentum().y(), + // (it + 1)->updatedState().globalMomentum().z()}; + // // calculate sigma_p on p = sqrt(p_x^2 + p_y^2 + p_z^2) as: + // // 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 + // for(int i = 3; i < 6; i++){ + // for(int j = 3; j < 6; j++){ + // // summing 2 * sigma_pij * p_i * p_j / p^2 + // // std::cout << "sigma_p_" << i << j << " = " << (it + 1)->updatedState().cartesianError().matrix()(i, j) << " " << "p_" << i << " = " << p[i-3] << " " << "p_" << j << " = " << p[j-3] << " "; + // sigma_p_cartesian += (it + 1)->updatedState().cartesianError().matrix()(i, j) * p[i-3] * p[j-3] / (it + 1)->updatedState().globalMomentum().mag2(); + // } + // } + + // sigma_p_cartesian = sqrt(sigma_p_cartesian); + + // sigma(p) using curvilinear error (on q/p) + float sigma_p = sqrt((it + 1)->updatedState().curvilinearError().matrix()(0, 0)) * (it + 1)->updatedState().globalMomentum().mag2(); + + trs.addSegment(layerpathlength, (it + 1)->updatedState().globalMomentum().mag2(), sigma_p); + LogTrace("TrackExtenderWithMTD") << "TSOS " << std::fixed << std::setw(4) << trs.size() << " R_i " << std::fixed << std::setw(14) << it->updatedState().globalPosition().perp() << " z_i " << std::fixed << std::setw(14) << it->updatedState().globalPosition().z() @@ -332,7 +399,14 @@ namespace { << std::setw(14) << (it + 1)->updatedState().globalPosition().z() << " p " << std::fixed << std::setw(14) << (it + 1)->updatedState().globalMomentum().mag() << " dp " << std::fixed << std::setw(14) - << (it + 1)->updatedState().globalMomentum().mag() - oldp; + << (it + 1)->updatedState().globalMomentum().mag() - oldp + << " sigma_p (from curvilinear) = " << std::fixed << std::setw(14) + << sigma_p + // << " sigma_p (cartesian) = " << std::fixed << std::setw(14) + // << sigma_p_cartesian + << " sigma_p/p = " << std::fixed << std::setw(14) + << sigma_p / (it + 1)->updatedState().globalMomentum().mag() * 100 << " %"; + oldp = (it + 1)->updatedState().globalMomentum().mag(); } @@ -345,12 +419,21 @@ namespace { validpropagation = false; } pathlength = pathlength1 + pathlength2; - trs.addSegment(pathlength2, tscblPCA.momentum().mag2()); + + float sigma_p = sqrt(tscblPCA.curvilinearError().matrix()(0, 0)) * tscblPCA.momentum().mag2(); + + trs.addSegment(pathlength2, tscblPCA.momentum().mag2(), sigma_p); + LogTrace("TrackExtenderWithMTD") << "TSOS " << std::fixed << std::setw(4) << trs.size() << " R_e " << std::fixed << std::setw(14) << tscblPCA.position().perp() << " z_e " << std::fixed << std::setw(14) << tscblPCA.position().z() << " p " << std::fixed << std::setw(14) << tscblPCA.momentum().mag() << " dp " << std::fixed << std::setw(14) - << tscblPCA.momentum().mag() - oldp; + << tscblPCA.momentum().mag() - oldp + << " sigma_p = " << std::fixed << std::setw(14) + << sigma_p + << " sigma_p/p = " << std::fixed << std::setw(14) + << sigma_p / tscblPCA.momentum().mag() * 100 << " %"; + return validpropagation; } @@ -459,7 +542,10 @@ class TrackExtenderWithMTDT : public edm::stream::EDProducer<> { float& sigmatmtdOut, float& tofpi, float& tofk, - float& tofp) const; + float& tofp, + float& sigmatofpi, + float& sigmatofk, + float& sigmatofp) const; reco::TrackExtra buildTrackExtra(const Trajectory& trajectory) const; string dumpLayer(const DetLayer* layer) const; @@ -481,6 +567,9 @@ class TrackExtenderWithMTDT : public edm::stream::EDProducer<> { edm::EDPutToken tofpiOrigTrkToken_; edm::EDPutToken tofkOrigTrkToken_; edm::EDPutToken tofpOrigTrkToken_; + edm::EDPutToken sigmatofpiOrigTrkToken_; + edm::EDPutToken sigmatofkOrigTrkToken_; + edm::EDPutToken sigmatofpOrigTrkToken_; edm::EDPutToken assocOrigTrkToken_; edm::EDGetTokenT tracksToken_; @@ -569,6 +658,9 @@ TrackExtenderWithMTDT::TrackExtenderWithMTDT(const ParameterSet tofpiOrigTrkToken_ = produces>("generalTrackTofPi"); tofkOrigTrkToken_ = produces>("generalTrackTofK"); tofpOrigTrkToken_ = produces>("generalTrackTofP"); + sigmatofpiOrigTrkToken_ = produces>("generalTrackSigmaTofPi"); + sigmatofkOrigTrkToken_ = produces>("generalTrackSigmaTofK"); + sigmatofpOrigTrkToken_ = produces>("generalTrackSigmaTofP"); assocOrigTrkToken_ = produces>("generalTrackassoc"); builderToken_ = esConsumes(edm::ESInputTag("", transientTrackBuilder_)); @@ -683,6 +775,9 @@ void TrackExtenderWithMTDT::produce(edm::Event& ev, const edm:: std::vector tofpiOrigTrkRaw; std::vector tofkOrigTrkRaw; std::vector tofpOrigTrkRaw; + std::vector sigmatofpiOrigTrkRaw; + std::vector sigmatofkOrigTrkRaw; + std::vector sigmatofpOrigTrkRaw; std::vector assocOrigTrkRaw; auto const tracksH = ev.getHandle(tracksToken_); @@ -727,6 +822,9 @@ void TrackExtenderWithMTDT::produce(edm::Event& ev, const edm:: LogTrace("TrackExtenderWithMTD") << "TrackExtenderWithMTD: extrapolating track " << itrack << " p/pT = " << track->p() << " " << track->pt() << " eta = " << track->eta(); + LogTrace("TrackExtenderWithMTD") << "TrackExtenderWithMTD: sigma_p = " + << sqrt(track->covariance()(0,0)) * track->p2() + << " sigma_p/p = " << sqrt(track->covariance()(0,0)) * track->p() * 100 << " %"; float trackVtxTime = 0.f; if (useVertex_) { @@ -803,12 +901,14 @@ void TrackExtenderWithMTDT::produce(edm::Event& ev, const edm:: const auto& trajwithmtd = mtdthits.empty() ? std::vector(1, trajs) : theTransformer->transform(ttrack, thits); float pMap = 0.f, betaMap = 0.f, t0Map = 0.f, sigmat0Map = -1.f, pathLengthMap = -1.f, tmtdMap = 0.f, - sigmatmtdMap = -1.f, tofpiMap = 0.f, tofkMap = 0.f, tofpMap = 0.f; + sigmatmtdMap = -1.f, tofpiMap = 0.f, tofkMap = 0.f, tofpMap = 0.f, sigmatofpiMap = -1.f, + sigmatofkMap = -1.f, sigmatofpMap = -1.f; int iMap = -1; for (const auto& trj : trajwithmtd) { const auto& thetrj = (updateTraj_ ? trj : trajs); - float pathLength = 0.f, tmtd = 0.f, sigmatmtd = -1.f, tofpi = 0.f, tofk = 0.f, tofp = 0.f; + float pathLength = 0.f, tmtd = 0.f, sigmatmtd = -1.f, tofpi = 0.f, tofk = 0.f, tofp = 0.f, sigmatofpi = -1.f, + sigmatofk = -1.f, sigmatofp = -1.f; LogTrace("TrackExtenderWithMTD") << "TrackExtenderWithMTD: refit track " << itrack << " p/pT = " << track->p() << " " << track->pt() << " eta = " << track->eta(); reco::Track result = buildTrack(track, @@ -823,7 +923,10 @@ void TrackExtenderWithMTDT::produce(edm::Event& ev, const edm:: sigmatmtd, tofpi, tofk, - tofp); + tofp, + sigmatofpi, + sigmatofk, + sigmatofp); if (result.ndof() >= 0) { /// setup the track extras reco::TrackExtra::TrajParams trajParams; @@ -856,6 +959,9 @@ void TrackExtenderWithMTDT::produce(edm::Event& ev, const edm:: tofpiMap = tofpi; tofkMap = tofk; tofpMap = tofp; + sigmatofpiMap = sigmatofpi; + sigmatofkMap = sigmatofk; + sigmatofpMap = sigmatofp; reco::TrackExtraRef extraRef(extrasRefProd, extras->size() - 1); backtrack.setExtra((updateExtra_ ? extraRef : track->extra())); for (unsigned ihit = hitsstart; ihit < hitsend; ++ihit) { @@ -864,8 +970,10 @@ void TrackExtenderWithMTDT::produce(edm::Event& ev, const edm:: npixBarrel.push_back(backtrack.hitPattern().numberOfValidPixelBarrelHits()); npixEndcap.push_back(backtrack.hitPattern().numberOfValidPixelEndcapHits()); LogTrace("TrackExtenderWithMTD") << "TrackExtenderWithMTD: tmtd " << tmtdMap << " +/- " << sigmatmtdMap - << " t0 " << t0Map << " +/- " << sigmat0Map << " tof pi/K/p " << tofpiMap - << " " << tofkMap << " " << tofpMap; + << " t0 " << t0Map << " +/- " << sigmat0Map << " tof pi/K/p " + << tofpiMap << "+/-" << fmt::format("{:0.2g}", sigmatofpiMap) << " (" << fmt::format("{:0.2g}", sigmatofpiMap/tofpiMap*100) << "%) " + << tofkMap << "+/-" << fmt::format("{:0.2g}", sigmatofkMap) << " (" << fmt::format("{:0.2g}", sigmatofkMap/tofkMap*100) << "%) " + << tofpMap << "+/-" << fmt::format("{:0.2g}", sigmatofpMap) << " (" << fmt::format("{:0.2g}", sigmatofpMap/tofpMap*100) << "%) "; } else { LogTrace("TrackExtenderWithMTD") << "Error in the MTD track refitting. This should not happen"; } @@ -881,6 +989,9 @@ void TrackExtenderWithMTDT::produce(edm::Event& ev, const edm:: tofpiOrigTrkRaw.push_back(tofpiMap); tofkOrigTrkRaw.push_back(tofkMap); tofpOrigTrkRaw.push_back(tofpMap); + sigmatofpiOrigTrkRaw.push_back(sigmatofpiMap); + sigmatofkOrigTrkRaw.push_back(sigmatofkMap); + sigmatofpOrigTrkRaw.push_back(sigmatofpMap); assocOrigTrkRaw.push_back(iMap); if (iMap == -1) { @@ -915,6 +1026,9 @@ void TrackExtenderWithMTDT::produce(edm::Event& ev, const edm:: fillValueMap(ev, tracksH, tofpiOrigTrkRaw, tofpiOrigTrkToken_); fillValueMap(ev, tracksH, tofkOrigTrkRaw, tofkOrigTrkToken_); fillValueMap(ev, tracksH, tofpOrigTrkRaw, tofpOrigTrkToken_); + fillValueMap(ev, tracksH, sigmatofpiOrigTrkRaw, sigmatofpiOrigTrkToken_); + fillValueMap(ev, tracksH, sigmatofkOrigTrkRaw, sigmatofkOrigTrkToken_); + fillValueMap(ev, tracksH, sigmatofpOrigTrkRaw, sigmatofpOrigTrkToken_); fillValueMap(ev, tracksH, assocOrigTrkRaw, assocOrigTrkToken_); } @@ -1176,7 +1290,11 @@ reco::Track TrackExtenderWithMTDT::buildTrack(const reco::Track float& sigmatmtdOut, float& tofpi, float& tofk, - float& tofp) const { + float& tofp, + float& sigmatofpi, + float& sigmatofk, + float& sigmatofp + ) const { TrajectoryStateClosestToBeamLine tscbl; bool tsbcl_status = getTrajectoryStateClosestToBeamLine(traj, bs, thePropagator, tscbl); @@ -1308,7 +1426,7 @@ reco::Track TrackExtenderWithMTDT::buildTrack(const reco::Track if (validmtd && validpropagation) { //here add the PID uncertainty for later use in the 1st step of 4D vtx reconstruction TrackTofPidInfo tofInfo = - computeTrackTofPidInfo(p.mag2(), pathlength, trs, thit, thiterror, 0.f, 0.f, true, TofCalc::kSegm); + computeTrackTofPidInfo(p.mag2(), pathlength, trs, thit, thiterror, 0.f, 0.f, true, TofCalc::kSegm, SigmaTofCalc::kCost); pathLengthOut = pathlength; // set path length if we've got a timing hit tmtdOut = thit; sigmatmtdOut = thiterror; @@ -1319,6 +1437,9 @@ reco::Track TrackExtenderWithMTDT::buildTrack(const reco::Track tofpi = tofInfo.dt_pi; tofk = tofInfo.dt_k; tofp = tofInfo.dt_p; + sigmatofpi = tofInfo.sigma_dt_pi; + sigmatofk = tofInfo.sigma_dt_k; + sigmatofp = tofInfo.sigma_dt_p; } } @@ -1426,4 +1547,4 @@ string TrackExtenderWithMTDT::dumpLayer(const DetLayer* layer) #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h" typedef TrackExtenderWithMTDT TrackExtenderWithMTD; -DEFINE_FWK_MODULE(TrackExtenderWithMTD); +DEFINE_FWK_MODULE(TrackExtenderWithMTD); \ No newline at end of file