Skip to content

Commit

Permalink
Fix in TPC MC cl.resoluion extraction
Browse files Browse the repository at this point in the history
  • Loading branch information
@shahor02
shahor02 committed Oct 28, 2024
1 parent f485cc0 commit 2c31473
Show file tree
Hide file tree
Showing 4 changed files with 138 additions and 52 deletions.
7 changes: 3 additions & 4 deletions Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TrackMCStudyConfig.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,10 +27,9 @@ struct TrackMCStudyConfig : o2::conf::ConfigurableParamHelper<TrackMCStudyConfig
float decayMotherMaxT = 1.0f; // max TOF in ns for mother particles to study
bool requireITSorTPCTrackRefs = true;
bool requireTopBottomRefs = false;
int minTPCRefsToExtractClRes = 4;
float rejectClustersResStat = 0.6;
float maxTPCRefExtrap = 2; // max dX to extrapolate the track ref when extrapolating track true posions
float maxTRefExtrapErr = 0.005;
int minTPCRefsToExtractClRes = 2;
float rejectClustersResStat = 0.;
float maxTPCRefExtrap = 2; // max dX to extrapolate the track ref when extrapolating track true posions
int decayPDG[5] = {310, 3122, 411, 421, -1}; // decays to study, must end by -1
O2ParamDef(TrackMCStudyConfig, "trmcconf");
};
Expand Down
78 changes: 61 additions & 17 deletions Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TrackMCStudyTypes.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include "CommonDataFormat/TimeStamp.h"
#include "ReconstructionDataFormats/PrimaryVertex.h"
#include <array>
#include <vector>

namespace o2::trackstudy
{
Expand Down Expand Up @@ -112,24 +113,20 @@ struct TrackFamily { // set of tracks related to the same MC label
ClassDefNV(TrackFamily, 1);
};

struct ClResTPC {
uint8_t sect = 0;
uint8_t row = 0;
uint8_t ncont = 0;
uint8_t flags = 0;
float snp = 0;
float tgl = 0;
float qmax = 0;
float qtot = 0;
float occ = 0;
std::array<float, 3> clPos{};
struct ClResTPCCont {
// contributor to TPC Cluster
std::array<float, 3> xyz{};
std::array<float, 3> below{};
std::array<float, 3> above{};
float snp = 0.;
float tgl = 0.;
float q2pt = 0.;
bool corrAttach = false;

int getNExt() const { return below[0] > 1. + above[0] > 1.; }
int getNExt() const { return (below[0] > 1.) + (above[0] > 1.); }

float getDY() const { return clPos[1] - getYRef(); }
float getDZ() const { return clPos[2] - getZRef(); }
float getDY() const { return xyz[1] - getYRef(); }
float getDZ() const { return xyz[2] - getZRef(); }

float getYRef() const
{
Expand Down Expand Up @@ -165,15 +162,62 @@ struct ClResTPC {
{
float adxA = 1e9, adxB = 1e9;
if (above[0] > 1.) {
adxA = clPos[0] - above[0];
adxA = xyz[0] - above[0];
}
if (below[0] > 1.) {
adxB = clPos[0] - below[0];
adxB = xyz[1] - below[0];
}
return std::abs(adxA) < std::abs(adxB) ? adxA : adxB;
}

ClassDefNV(ClResTPC, 1);
float getDXMax() const
{
float adxA = 0, adxB = 0;
if (above[0] > 1.) {
adxA = xyz[0] - above[0];
}
if (below[0] > 1.) {
adxB = xyz[0] - below[0];
}
return std::abs(adxA) > std::abs(adxB) ? adxA : adxB;
}

float getEY() const { return getNExt() > 1 ? below[1] - above[1] : -999; }
float getEZ() const { return getNExt() > 1 ? below[2] - above[2] : -999; }

ClassDefNV(ClResTPCCont, 1);
};

struct ClResTPC {
uint8_t sect = 0;
uint8_t row = 0;
uint8_t ncont = 0;
uint8_t flags = 0;
float qmax = 0;
float qtot = 0;
float occ = 0;

std::vector<ClResTPCCont> contTracks;
int getNCont() const { return contTracks.size(); }

float getDY(int i) const { return i < getNCont() ? contTracks[i].getDY() : -999.; }
float getDZ(int i) const { return i < getNCont() ? contTracks[i].getDZ() : -999.; }
float getYRef(int i) const { return i < getNCont() ? contTracks[i].getYRef() : -999.; }
float getZRef(int i) const { return i < getNCont() ? contTracks[i].getZRef() : -999.; }
float getDXMin(int i) const { return i < getNCont() ? contTracks[i].getDXMin() : -999.; }
float getDXMax(int i) const { return i < getNCont() ? contTracks[i].getDXMax() : -999.; }
float getEY(int i) const { return i < getNCont() ? contTracks[i].getEY() : -999.; }
float getEZ(int i) const { return i < getNCont() ? contTracks[i].getEZ() : -999.; }

void sortCont()
{
std::sort(contTracks.begin(), contTracks.end(), [](const ClResTPCCont& a, const ClResTPCCont& b) {
float dya = a.getDY(), dyb = b.getDY(), dza = a.getDZ(), dzb = b.getDZ();
return dya * dya + dza * dza < dyb * dyb + dzb * dzb;
});
}

ClassDefNV(ClResTPC, 2);
};

struct RecPV {
Expand Down
2 changes: 2 additions & 0 deletions Detectors/GlobalTrackingWorkflow/study/src/GlobalTrackingStudyLinkDef.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,5 +34,7 @@
#pragma link C++ class o2::trackstudy::MCVertex + ;
#pragma link C++ class std::vector < o2::trackstudy::MCVertex> + ;
#pragma link C++ class o2::trackstudy::ClResTPC + ;
#pragma link C++ class o2::trackstudy::ClResTPCCont + ;
#pragma link C++ class std::vector < o2::trackstudy::ClResTPCCont> + ;

#endif
103 changes: 72 additions & 31 deletions Detectors/GlobalTrackingWorkflow/study/src/TrackMCStudy.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -263,6 +263,26 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
return -1;
};

auto flagTPCClusters = [&recoData](const o2::tpc::TrackTPC& trc, o2::MCCompLabel lbTrc) {
if (recoData.inputsTPCclusters) {
const auto clRefs = recoData.getTPCTracksClusterRefs();
const auto* TPCClMClab = recoData.inputsTPCclusters->clusterIndex.clustersMCTruth;
const auto& TPCClusterIdxStruct = recoData.inputsTPCclusters->clusterIndex;
for (int ic = 0; ic < trc.getNClusterReferences(); ic++) {
uint8_t clSect = 0, clRow = 0;
uint32_t clIdx = 0;
trc.getClusterReference(clRefs, ic, clSect, clRow, clIdx);
auto labels = TPCClMClab->getLabels(clIdx + TPCClusterIdxStruct.clusterOffset[clSect][clRow]);
for (auto& lbl : labels) {
if (lbl == lbTrc) {
const_cast<o2::MCCompLabel&>(lbl).setFakeFlag(true); // actually, in this way we are flagging that this cluster was correctly attached
break;
}
}
}
}
};

{
const auto* digconst = mcReader.getDigitizationContext();
const auto& mcEvRecords = digconst->getEventRecords(false);
Expand Down Expand Up @@ -421,13 +441,11 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
}
}

// collect ITS/TPC cluster info for selected MC particles
if (params.minTPCRefsToExtractClRes > 0) {
LOGP(info, "collected {} MC tracks", mSelMCTracks.size());
if (params.minTPCRefsToExtractClRes > 0) { // prepare MC trackrefs for TPC
processTPCTrackRefs();
}
fillMCClusterInfo(recoData);

LOGP(info, "collected {} MC tracks", mSelMCTracks.size());
int mcnt = 0;
for (auto& entry : mSelMCTracks) {
auto& trackFam = entry.second;
Expand Down Expand Up @@ -489,6 +507,7 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
const auto& trtpc = recoData.getTPCTrack(gidSet[GTrackID::TPC]);
tref.nClTPC = trtpc.getNClusters();
tref.lowestPadRow = getLowestPadrow(trtpc);
flagTPCClusters(trtpc, entry.first);
if (trackFam.entTPC < 0) {
trackFam.entTPC = tcnt;
trackFam.tpcT0 = trtpc.getTime0();
Expand Down Expand Up @@ -533,6 +552,10 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
auto& trackFam = entry.second;
(*mDBGOut) << "tracks" << "tr=" << trackFam << "\n";
}

// collect ITS/TPC cluster info for selected MC particles
fillMCClusterInfo(recoData);

// decays
std::vector<TrackFamily> decFam;
for (int id = 0; id < mNCheckDecays; id++) {
Expand Down Expand Up @@ -620,6 +643,8 @@ void TrackMCStudy::fillMCClusterInfo(const o2::globaltracking::RecoContainer& re
const auto& TPCClusterIdxStruct = recoData.inputsTPCclusters->clusterIndex;
const auto* TPCClMClab = recoData.inputsTPCclusters->clusterIndex.clustersMCTruth;
const auto& params = o2::trackstudy::TrackMCStudyConfig::Instance();

ClResTPC clRes{};
for (uint8_t sector = 0; sector < 36; sector++) {
for (uint8_t row = 0; row < 152; row++) {
unsigned int offs = TPCClusterIdxStruct.clusterOffset[sector][row];
Expand All @@ -632,7 +657,12 @@ void TrackMCStudy::fillMCClusterInfo(const o2::globaltracking::RecoContainer& re
}
ncontLb++;
}
for (const auto& lbl : labels) {
const auto& clus = TPCClusterIdxStruct.clusters[sector][row][icl0];
clRes.contTracks.clear();
bool doClusRes = (params.minTPCRefsToExtractClRes > 0) && (params.rejectClustersResStat <= 0. || gRandom->Rndm() < params.rejectClustersResStat);
for (auto lbl : labels) {
bool corrAttach = lbl.isFake(); // was this flagged in the flagTPCClusters called from process ?
lbl.setFakeFlag(false);
auto entry = mSelMCTracks.find(lbl);
if (entry == mSelMCTracks.end()) { // not selected
continue;
Expand All @@ -654,29 +684,29 @@ void TrackMCStudy::fillMCClusterInfo(const o2::globaltracking::RecoContainer& re
mctr.nTPCClShared++;
}
// try to extract ideal track position
if (params.minTPCRefsToExtractClRes > 0) {
if (doClusRes) {
auto entTRefIDsIt = mSelTRefIdx.find(lbl);
if (entTRefIDsIt == mSelTRefIdx.end()) {
continue;
}
float xc, yc, zc;
mTPCCorrMapsLoader.Transform(sector, row, clus.getPad(), clus.getTime(), xc, yc, zc, mctr.bcInTF / 8.); // nominal time of the track

const auto& entTRefIDs = entTRefIDsIt->second;
// find bracketing TRef params
int entIDBelow = -1, entIDAbove = -1;
float xBelow = -1e6, xAbove = 1e6;
const auto& clus = TPCClusterIdxStruct.clusters[sector][row][icl0];
float xc, yc, zc;
mTPCCorrMapsLoader.Transform(sector, row, clus.getPad(), clus.getTime(), xc, yc, zc, mctr.bcInTF / 8.); // nominal time of the track

for (int entID = entTRefIDs.first; entID < entTRefIDs.second; entID++) {
const auto& refTr = mSelTRefs[entID];
if (refTr.getUserField() != sector % 18) {
continue;
}
if (refTr.getX() < xc && refTr.getX() > xBelow && refTr.getX() > xc - params.maxTPCRefExtrap) {
if ((refTr.getX() < xc) && (refTr.getX() > xBelow) && (refTr.getX() > xc - params.maxTPCRefExtrap)) {
xBelow = refTr.getX();
entIDBelow = entID;
}
if (refTr.getX() > xc && refTr.getX() < xAbove && refTr.getX() < xc + params.maxTPCRefExtrap) {
if ((refTr.getX() > xc) && (refTr.getX() < xAbove) && (refTr.getX() < xc + params.maxTPCRefExtrap)) {
xAbove = refTr.getX();
entIDAbove = entID;
}
Expand All @@ -688,42 +718,53 @@ void TrackMCStudy::fillMCClusterInfo(const o2::globaltracking::RecoContainer& re
o2::track::TrackPar tparAbove, tparBelow;
bool okBelow = entIDBelow >= 0 && prop->PropagateToXBxByBz((tparBelow = mSelTRefs[entIDBelow]), xc, 0.99, 2.);
bool okAbove = entIDAbove >= 0 && prop->PropagateToXBxByBz((tparAbove = mSelTRefs[entIDAbove]), xc, 0.99, 2.);
if ((!okBelow && !okAbove) || (params.requireTopBottomRefs && (!okBelow || !okAbove)) || (params.rejectClustersResStat > 0. && gRandom->Rndm() < params.rejectClustersResStat)) {
if ((!okBelow && !okAbove) || (params.requireTopBottomRefs && (!okBelow || !okAbove))) {
continue;
}
ClResTPC clRes{};

int nmeas = 0;
auto& clCont = clRes.contTracks.emplace_back();
clCont.corrAttach = corrAttach;
if (okBelow) {
clRes.below = {mSelTRefs[entIDBelow].getX(), tparBelow.getY(), tparBelow.getZ()};
clRes.snp += tparBelow.getSnp();
clRes.tgl += tparBelow.getTgl();
clCont.below = {mSelTRefs[entIDBelow].getX(), tparBelow.getY(), tparBelow.getZ()};
clCont.snp += tparBelow.getSnp();
clCont.tgl += tparBelow.getTgl();
clCont.q2pt += tparBelow.getQ2Pt();
nmeas++;
}
if (okAbove) {
clRes.above = {mSelTRefs[entIDAbove].getX(), tparAbove.getY(), tparAbove.getZ()};
clRes.snp += tparAbove.getSnp();
clRes.tgl += tparBelow.getTgl();
clCont.above = {mSelTRefs[entIDAbove].getX(), tparAbove.getY(), tparAbove.getZ()};
clCont.snp += tparAbove.getSnp();
clCont.tgl += tparAbove.getTgl();
clCont.q2pt += tparAbove.getQ2Pt();
nmeas++;
}
if (nmeas) {
if (clRes.contTracks.size() == 1) {
int occBin = mctr.bcInTF / 8 * mNTPCOccBinLengthInv;
clRes.occ = occBin < 0 ? mTBinClOcc[0] : (occBin >= mTBinClOcc.size() ? mTBinClOcc.back() : mTBinClOcc[occBin]);
}
clCont.xyz = {xc, yc, zc};
if (nmeas > 1) {
clRes.snp *= 0.5;
clRes.tgl *= 0.5;
clCont.snp *= 0.5;
clCont.tgl *= 0.5;
clCont.q2pt *= 0.5;
}
clRes.clPos = {xc, yc, zc};
clRes.sect = sector;
clRes.row = row;
clRes.qtot = clus.getQtot();
clRes.qmax = clus.getQmax();
clRes.flags = clus.getFlags();
clRes.ncont = ncontLb;
int occBin = mctr.bcInTF / 8 * mNTPCOccBinLengthInv;
clRes.occ = occBin < 0 ? mTBinClOcc[0] : (occBin >= mTBinClOcc.size() ? mTBinClOcc.back() : mTBinClOcc[occBin]);
(*mDBGOut) << "clres" << "clr=" << clRes << "\n";
} else {
clRes.contTracks.pop_back();
}
}
}
if (clRes.getNCont()) {
clRes.sect = sector;
clRes.row = row;
clRes.qtot = clus.getQtot();
clRes.qmax = clus.getQmax();
clRes.flags = clus.getFlags();
clRes.ncont = ncontLb;
clRes.sortCont();
(*mDBGOut) << "clres" << "clr=" << clRes << "\n";
}
}
}
}
Expand Down

0 comments on commit 2c31473

Please sign in to comment.