Skip to content

Commit

Permalink
Squashed commit of the following:
Browse files Browse the repository at this point in the history 
commit a7454c6
Author: noepalm <palmeri.1871005@studenti.uniroma1.it>
Date:   Mon Feb 12 14:22:51 2024 +0100

    Updated binning for sigma(TOF) validation histograms

    Updated binning for sigma(TOF) validation histograms

    Updated binning for sigma(TOF) validation histograms

commit e0dc3be
Author: noepalm <palmeri.1871005@studenti.uniroma1.it>
Date:   Wed Feb 7 16:40:42 2024 +0100

    Code checks && format (#2)

commit 7c70ed8
Author: noepalm <palmeri.1871005@studenti.uniroma1.it>
Date:   Wed Feb 7 15:57:59 2024 +0100

    Removed debug ifdef condition on print to match previous implementation

commit 2e40a9f
Author: noepalm <palmeri.1871005@studenti.uniroma1.it>
Date:   Wed Feb 7 15:54:41 2024 +0100

    Removed useless comments, files + fixed input file path for standalone MTD validation cfg file

commit 199ad72
Author: noepalm <palmeri.1871005@studenti.uniroma1.it>
Date:   Wed Feb 7 14:41:57 2024 +0100

    Code checks && code format

commit 6ec3c46
Author: noepalm <palmeri.1871005@studenti.uniroma1.it>
Date:   Tue Feb 6 09:56:49 2024 +0100

    Changed SigmaTof fill location in validation

commit 6a0525e
Author: noepalm <palmeri.1871005@studenti.uniroma1.it>
Date:   Mon Feb 5 18:11:27 2024 +0100

    Change input file to validation

commit cb7db5b
Author: noepalm <palmeri.1871005@studenti.uniroma1.it>
Date:   Mon Feb 5 18:07:19 2024 +0100

    Saving sigmaTOF in ps

commit 00230ac
Author: noepalm <palmeri.1871005@studenti.uniroma1.it>
Date:   Mon Feb 5 15:32:11 2024 +0100

    Removed useless prints + created member vector in TrackSegments to allocate memory for sigmaTof once

commit fdd6e31
Author: noepalm <palmeri.1871005@studenti.uniroma1.it>
Date:   Tue Jan 9 16:08:46 2024 +0100

    Added sigma(p) to TrackSegments; added sigma(TOF) computation, now saved to event

commit 7d488e0
Author: noepalm <palmeri.1871005@studenti.uniroma1.it>
Date:   Wed Dec 6 14:24:26 2023 +0100

    testing files

commit 88f9847
Author: noepalm <palmeri.1871005@studenti.uniroma1.it>
Date:   Mon Feb 12 12:42:15 2024 +0100

    Replaced isnan, isinf with edm::isNotFinite in MTD vertex validation
  • Loading branch information
@noepalm
noepalm committed Feb 12, 2024
1 parent 7d5cb36 commit b2a7305
Show file tree
Hide file tree
Showing 3 changed files with 197 additions and 19 deletions.
155 changes: 139 additions & 16 deletions RecoMTD/TrackExtender/plugins/TrackExtenderWithMTD.cc
View file
Expand Up @@ -96,15 +96,19 @@ namespace {

class TrackSegments {
public:
TrackSegments() = default;
TrackSegments() {
sigmaTofs_.reserve(30); // observed upper limit on nSegments
};

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_;
}
Expand All @@ -118,11 +122,49 @@ namespace {

LogTrace("TrackExtenderWithMTD") << " TOF Segment # " << iSeg + 1 << " p = " << std::sqrt(segmentMom2_[iSeg])
<< " tof = " << tof;

#ifdef EDM_ML_DEBUG
float sigma_tof = segmentPathOvc_[iSeg] * segmentSigmaMom_[iSeg] /
(segmentMom2_[iSeg] * sqrt(segmentMom2_[iSeg] + 1 / mass_inv2) * mass_inv2);

LogTrace("TrackExtenderWithMTD") << "TOF Segment # " << iSeg + 1 << std::fixed << std::setw(6)
<< " tof segment = " << segmentPathOvc_[iSeg] / beta << std::scientific
<< "+/- " << sigma_tof << std::fixed
<< "(rel. err. = " << sigma_tof / (segmentPathOvc_[iSeg] / beta) * 100
<< " %)";
#endif
}

return tof;
}

inline float computeSigmaTof(float mass_inv2) {
float sigmatof = 0.;

// remove previously calculated sigmaTofs
sigmaTofs_.clear();

// compute sigma(tof) on each segment first by propagating sigma(p)
// also add diagonal terms to sigmatof
float sigma = 0.;
for (uint32_t iSeg = 0; iSeg < nSegment_; iSeg++) {
sigma = segmentPathOvc_[iSeg] * segmentSigmaMom_[iSeg] /
(segmentMom2_[iSeg] * sqrt(segmentMom2_[iSeg] + 1 / mass_inv2) * mass_inv2);
sigmaTofs_.push_back(sigma);

sigmatof += sigma * sigma;
}

// compute sigma on sum of tofs assuming full correlation between segments
for (uint32_t iSeg = 0; iSeg < nSegment_; iSeg++) {
for (uint32_t jSeg = iSeg + 1; jSeg < nSegment_; jSeg++) {
sigmatof += 2 * sigmaTofs_[iSeg] * sigmaTofs_[jSeg];
}
}

return sqrt(sigmatof);
}

inline uint32_t size() const { return nSegment_; }

inline uint32_t removeFirstSegment() {
Expand All @@ -144,6 +186,9 @@ namespace {
uint32_t nSegment_ = 0;
std::vector<float> segmentPathOvc_;
std::vector<float> segmentMom2_;
std::vector<float> segmentSigmaMom_;

std::vector<float> sigmaTofs_;
};

struct TrackTofPidInfo {
Expand All @@ -164,21 +209,25 @@ 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;
float prob_p;
};

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

const TrackTofPidInfo computeTrackTofPidInfo(float magp2,
float length,
Expand All @@ -188,7 +237,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;
Expand Down Expand Up @@ -218,17 +268,36 @@ 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;
};

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);

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);

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);
Expand Down Expand Up @@ -323,7 +392,13 @@ namespace {
validpropagation = false;
}
pathlength1 += layerpathlength;
trs.addSegment(layerpathlength, (it + 1)->updatedState().globalMomentum().mag2());

// sigma(p) from 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()
Expand All @@ -345,12 +420,19 @@ 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;
}

Expand Down Expand Up @@ -459,7 +541,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;
Expand All @@ -481,6 +566,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<InputCollection> tracksToken_;
Expand Down Expand Up @@ -569,6 +657,9 @@ TrackExtenderWithMTDT<TrackCollection>::TrackExtenderWithMTDT(const ParameterSet
tofpiOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackTofPi");
tofkOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackTofK");
tofpOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackTofP");
sigmatofpiOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackSigmaTofPi");
sigmatofkOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackSigmaTofK");
sigmatofpOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackSigmaTofP");
assocOrigTrkToken_ = produces<edm::ValueMap<int>>("generalTrackassoc");

builderToken_ = esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", transientTrackBuilder_));
Expand Down Expand Up @@ -683,6 +774,9 @@ void TrackExtenderWithMTDT<TrackCollection>::produce(edm::Event& ev, const edm::
std::vector<float> tofpiOrigTrkRaw;
std::vector<float> tofkOrigTrkRaw;
std::vector<float> tofpOrigTrkRaw;
std::vector<float> sigmatofpiOrigTrkRaw;
std::vector<float> sigmatofkOrigTrkRaw;
std::vector<float> sigmatofpOrigTrkRaw;
std::vector<int> assocOrigTrkRaw;

auto const tracksH = ev.getHandle(tracksToken_);
Expand Down Expand Up @@ -727,6 +821,9 @@ void TrackExtenderWithMTDT<TrackCollection>::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_) {
Expand Down Expand Up @@ -803,12 +900,14 @@ void TrackExtenderWithMTDT<TrackCollection>::produce(edm::Event& ev, const edm::
const auto& trajwithmtd =
mtdthits.empty() ? std::vector<Trajectory>(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,
Expand All @@ -823,7 +922,10 @@ void TrackExtenderWithMTDT<TrackCollection>::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;
Expand Down Expand Up @@ -856,6 +958,9 @@ void TrackExtenderWithMTDT<TrackCollection>::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) {
Expand All @@ -865,7 +970,12 @@ void TrackExtenderWithMTDT<TrackCollection>::produce(edm::Event& ev, const edm::
npixEndcap.push_back(backtrack.hitPattern().numberOfValidPixelEndcapHits());
LogTrace("TrackExtenderWithMTD") << "TrackExtenderWithMTD: tmtd " << tmtdMap << " +/- " << sigmatmtdMap
<< " t0 " << t0Map << " +/- " << sigmat0Map << " tof pi/K/p " << tofpiMap
<< " " << tofkMap << " " << tofpMap;
<< "+/-" << 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";
}
Expand All @@ -881,6 +991,9 @@ void TrackExtenderWithMTDT<TrackCollection>::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) {
Expand Down Expand Up @@ -915,6 +1028,9 @@ void TrackExtenderWithMTDT<TrackCollection>::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_);
}

Expand Down Expand Up @@ -1176,7 +1292,10 @@ reco::Track TrackExtenderWithMTDT<TrackCollection>::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);

Expand Down Expand Up @@ -1307,8 +1426,9 @@ reco::Track TrackExtenderWithMTDT<TrackCollection>::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);
TrackTofPidInfo tofInfo = 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;
Expand All @@ -1319,6 +1439,9 @@ reco::Track TrackExtenderWithMTDT<TrackCollection>::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;
}
}

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

DEFINE_FWK_MODULE(TrackExtenderWithMTD);
DEFINE_FWK_MODULE(TrackExtenderWithMTD);

0 comments on commit b2a7305

Please sign in to comment.