Skip to content

Commit

Permalink
Added sigma(p) to TrackSegments; added sigma(TOF) computation, now sa…
Browse files Browse the repository at this point in the history
…ved to event
  • Loading branch information
noepalm committed Feb 12, 2024
1 parent 7d488e0 commit fdd6e31
Showing 1 changed file with 139 additions and 18 deletions.
157 changes: 139 additions & 18 deletions RecoMTD/TrackExtender/plugins/TrackExtenderWithMTD.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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_;
}
Expand All @@ -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_; }

Expand All @@ -144,6 +163,7 @@ namespace {
uint32_t nSegment_ = 0;
std::vector<float> segmentPathOvc_;
std::vector<float> segmentMom2_;
std::vector<float> segmentSigmaMom_;
};

struct TrackTofPidInfo {
Expand All @@ -164,21 +184,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 +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;
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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()
Expand All @@ -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();
}

Expand All @@ -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;
}

Expand Down Expand Up @@ -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;
Expand All @@ -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<InputCollection> tracksToken_;
Expand Down Expand Up @@ -569,6 +658,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 +775,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 +822,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 +901,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 +923,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 +959,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 @@ -864,8 +970,10 @@ void TrackExtenderWithMTDT<TrackCollection>::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";
}
Expand All @@ -881,6 +989,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 +1026,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 +1290,11 @@ 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 @@ -1308,7 +1426,7 @@ 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);
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 +1437,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 +1547,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 fdd6e31

Please sign in to comment.