Skip to content

Commit 2c31473

Browse files
shahor02shahor02
committed
Fix in TPC MC cl.resoluion extraction
1 parent f485cc0 commit 2c31473

File tree

4 files changed

+138
-52
lines changed

4 files changed

+138
-52
lines changed

‎Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TrackMCStudyConfig.h‎

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -27,10 +27,9 @@ struct TrackMCStudyConfig : o2::conf::ConfigurableParamHelper<TrackMCStudyConfig
2727
float decayMotherMaxT = 1.0f; // max TOF in ns for mother particles to study
2828
bool requireITSorTPCTrackRefs = true;
2929
bool requireTopBottomRefs = false;
30-
int minTPCRefsToExtractClRes = 4;
31-
float rejectClustersResStat = 0.6;
32-
float maxTPCRefExtrap = 2; // max dX to extrapolate the track ref when extrapolating track true posions
33-
float maxTRefExtrapErr = 0.005;
30+
int minTPCRefsToExtractClRes = 2;
31+
float rejectClustersResStat = 0.;
32+
float maxTPCRefExtrap = 2; // max dX to extrapolate the track ref when extrapolating track true posions
3433
int decayPDG[5] = {310, 3122, 411, 421, -1}; // decays to study, must end by -1
3534
O2ParamDef(TrackMCStudyConfig, "trmcconf");
3635
};

‎Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TrackMCStudyTypes.h‎

Lines changed: 61 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@
2020
#include "CommonDataFormat/TimeStamp.h"
2121
#include "ReconstructionDataFormats/PrimaryVertex.h"
2222
#include <array>
23+
#include <vector>
2324

2425
namespace o2::trackstudy
2526
{
@@ -112,24 +113,20 @@ struct TrackFamily { // set of tracks related to the same MC label
112113
ClassDefNV(TrackFamily, 1);
113114
};
114115

115-
struct ClResTPC {
116-
uint8_t sect = 0;
117-
uint8_t row = 0;
118-
uint8_t ncont = 0;
119-
uint8_t flags = 0;
120-
float snp = 0;
121-
float tgl = 0;
122-
float qmax = 0;
123-
float qtot = 0;
124-
float occ = 0;
125-
std::array<float, 3> clPos{};
116+
struct ClResTPCCont {
117+
// contributor to TPC Cluster
118+
std::array<float, 3> xyz{};
126119
std::array<float, 3> below{};
127120
std::array<float, 3> above{};
121+
float snp = 0.;
122+
float tgl = 0.;
123+
float q2pt = 0.;
124+
bool corrAttach = false;
128125

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

131-
float getDY() const { return clPos[1] - getYRef(); }
132-
float getDZ() const { return clPos[2] - getZRef(); }
128+
float getDY() const { return xyz[1] - getYRef(); }
129+
float getDZ() const { return xyz[2] - getZRef(); }
133130

134131
float getYRef() const
135132
{
@@ -165,15 +162,62 @@ struct ClResTPC {
165162
{
166163
float adxA = 1e9, adxB = 1e9;
167164
if (above[0] > 1.) {
168-
adxA = clPos[0] - above[0];
165+
adxA = xyz[0] - above[0];
169166
}
170167
if (below[0] > 1.) {
171-
adxB = clPos[0] - below[0];
168+
adxB = xyz[1] - below[0];
172169
}
173170
return std::abs(adxA) < std::abs(adxB) ? adxA : adxB;
174171
}
175172

176-
ClassDefNV(ClResTPC, 1);
173+
float getDXMax() const
174+
{
175+
float adxA = 0, adxB = 0;
176+
if (above[0] > 1.) {
177+
adxA = xyz[0] - above[0];
178+
}
179+
if (below[0] > 1.) {
180+
adxB = xyz[0] - below[0];
181+
}
182+
return std::abs(adxA) > std::abs(adxB) ? adxA : adxB;
183+
}
184+
185+
float getEY() const { return getNExt() > 1 ? below[1] - above[1] : -999; }
186+
float getEZ() const { return getNExt() > 1 ? below[2] - above[2] : -999; }
187+
188+
ClassDefNV(ClResTPCCont, 1);
189+
};
190+
191+
struct ClResTPC {
192+
uint8_t sect = 0;
193+
uint8_t row = 0;
194+
uint8_t ncont = 0;
195+
uint8_t flags = 0;
196+
float qmax = 0;
197+
float qtot = 0;
198+
float occ = 0;
199+
200+
std::vector<ClResTPCCont> contTracks;
201+
int getNCont() const { return contTracks.size(); }
202+
203+
float getDY(int i) const { return i < getNCont() ? contTracks[i].getDY() : -999.; }
204+
float getDZ(int i) const { return i < getNCont() ? contTracks[i].getDZ() : -999.; }
205+
float getYRef(int i) const { return i < getNCont() ? contTracks[i].getYRef() : -999.; }
206+
float getZRef(int i) const { return i < getNCont() ? contTracks[i].getZRef() : -999.; }
207+
float getDXMin(int i) const { return i < getNCont() ? contTracks[i].getDXMin() : -999.; }
208+
float getDXMax(int i) const { return i < getNCont() ? contTracks[i].getDXMax() : -999.; }
209+
float getEY(int i) const { return i < getNCont() ? contTracks[i].getEY() : -999.; }
210+
float getEZ(int i) const { return i < getNCont() ? contTracks[i].getEZ() : -999.; }
211+
212+
void sortCont()
213+
{
214+
std::sort(contTracks.begin(), contTracks.end(), [](const ClResTPCCont& a, const ClResTPCCont& b) {
215+
float dya = a.getDY(), dyb = b.getDY(), dza = a.getDZ(), dzb = b.getDZ();
216+
return dya * dya + dza * dza < dyb * dyb + dzb * dzb;
217+
});
218+
}
219+
220+
ClassDefNV(ClResTPC, 2);
177221
};
178222

179223
struct RecPV {

‎Detectors/GlobalTrackingWorkflow/study/src/GlobalTrackingStudyLinkDef.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -34,5 +34,7 @@
3434
#pragma link C++ class o2::trackstudy::MCVertex + ;
3535
#pragma link C++ class std::vector < o2::trackstudy::MCVertex> + ;
3636
#pragma link C++ class o2::trackstudy::ClResTPC + ;
37+
#pragma link C++ class o2::trackstudy::ClResTPCCont + ;
38+
#pragma link C++ class std::vector < o2::trackstudy::ClResTPCCont> + ;
3739

3840
#endif

‎Detectors/GlobalTrackingWorkflow/study/src/TrackMCStudy.cxx‎

Lines changed: 72 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -263,6 +263,26 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
263263
return -1;
264264
};
265265

266+
auto flagTPCClusters = [&recoData](const o2::tpc::TrackTPC& trc, o2::MCCompLabel lbTrc) {
267+
if (recoData.inputsTPCclusters) {
268+
const auto clRefs = recoData.getTPCTracksClusterRefs();
269+
const auto* TPCClMClab = recoData.inputsTPCclusters->clusterIndex.clustersMCTruth;
270+
const auto& TPCClusterIdxStruct = recoData.inputsTPCclusters->clusterIndex;
271+
for (int ic = 0; ic < trc.getNClusterReferences(); ic++) {
272+
uint8_t clSect = 0, clRow = 0;
273+
uint32_t clIdx = 0;
274+
trc.getClusterReference(clRefs, ic, clSect, clRow, clIdx);
275+
auto labels = TPCClMClab->getLabels(clIdx + TPCClusterIdxStruct.clusterOffset[clSect][clRow]);
276+
for (auto& lbl : labels) {
277+
if (lbl == lbTrc) {
278+
const_cast<o2::MCCompLabel&>(lbl).setFakeFlag(true); // actually, in this way we are flagging that this cluster was correctly attached
279+
break;
280+
}
281+
}
282+
}
283+
}
284+
};
285+
266286
{
267287
const auto* digconst = mcReader.getDigitizationContext();
268288
const auto& mcEvRecords = digconst->getEventRecords(false);
@@ -421,13 +441,11 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
421441
}
422442
}
423443

424-
// collect ITS/TPC cluster info for selected MC particles
425-
if (params.minTPCRefsToExtractClRes > 0) {
444+
LOGP(info, "collected {} MC tracks", mSelMCTracks.size());
445+
if (params.minTPCRefsToExtractClRes > 0) { // prepare MC trackrefs for TPC
426446
processTPCTrackRefs();
427447
}
428-
fillMCClusterInfo(recoData);
429448

430-
LOGP(info, "collected {} MC tracks", mSelMCTracks.size());
431449
int mcnt = 0;
432450
for (auto& entry : mSelMCTracks) {
433451
auto& trackFam = entry.second;
@@ -489,6 +507,7 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
489507
const auto& trtpc = recoData.getTPCTrack(gidSet[GTrackID::TPC]);
490508
tref.nClTPC = trtpc.getNClusters();
491509
tref.lowestPadRow = getLowestPadrow(trtpc);
510+
flagTPCClusters(trtpc, entry.first);
492511
if (trackFam.entTPC < 0) {
493512
trackFam.entTPC = tcnt;
494513
trackFam.tpcT0 = trtpc.getTime0();
@@ -533,6 +552,10 @@ void TrackMCStudy::process(const o2::globaltracking::RecoContainer& recoData)
533552
auto& trackFam = entry.second;
534553
(*mDBGOut) << "tracks" << "tr=" << trackFam << "\n";
535554
}
555+
556+
// collect ITS/TPC cluster info for selected MC particles
557+
fillMCClusterInfo(recoData);
558+
536559
// decays
537560
std::vector<TrackFamily> decFam;
538561
for (int id = 0; id < mNCheckDecays; id++) {
@@ -620,6 +643,8 @@ void TrackMCStudy::fillMCClusterInfo(const o2::globaltracking::RecoContainer& re
620643
const auto& TPCClusterIdxStruct = recoData.inputsTPCclusters->clusterIndex;
621644
const auto* TPCClMClab = recoData.inputsTPCclusters->clusterIndex.clustersMCTruth;
622645
const auto& params = o2::trackstudy::TrackMCStudyConfig::Instance();
646+
647+
ClResTPC clRes{};
623648
for (uint8_t sector = 0; sector < 36; sector++) {
624649
for (uint8_t row = 0; row < 152; row++) {
625650
unsigned int offs = TPCClusterIdxStruct.clusterOffset[sector][row];
@@ -632,7 +657,12 @@ void TrackMCStudy::fillMCClusterInfo(const o2::globaltracking::RecoContainer& re
632657
}
633658
ncontLb++;
634659
}
635-
for (const auto& lbl : labels) {
660+
const auto& clus = TPCClusterIdxStruct.clusters[sector][row][icl0];
661+
clRes.contTracks.clear();
662+
bool doClusRes = (params.minTPCRefsToExtractClRes > 0) && (params.rejectClustersResStat <= 0. || gRandom->Rndm() < params.rejectClustersResStat);
663+
for (auto lbl : labels) {
664+
bool corrAttach = lbl.isFake(); // was this flagged in the flagTPCClusters called from process ?
665+
lbl.setFakeFlag(false);
636666
auto entry = mSelMCTracks.find(lbl);
637667
if (entry == mSelMCTracks.end()) { // not selected
638668
continue;
@@ -654,29 +684,29 @@ void TrackMCStudy::fillMCClusterInfo(const o2::globaltracking::RecoContainer& re
654684
mctr.nTPCClShared++;
655685
}
656686
// try to extract ideal track position
657-
if (params.minTPCRefsToExtractClRes > 0) {
687+
if (doClusRes) {
658688
auto entTRefIDsIt = mSelTRefIdx.find(lbl);
659689
if (entTRefIDsIt == mSelTRefIdx.end()) {
660690
continue;
661691
}
692+
float xc, yc, zc;
693+
mTPCCorrMapsLoader.Transform(sector, row, clus.getPad(), clus.getTime(), xc, yc, zc, mctr.bcInTF / 8.); // nominal time of the track
694+
662695
const auto& entTRefIDs = entTRefIDsIt->second;
663696
// find bracketing TRef params
664697
int entIDBelow = -1, entIDAbove = -1;
665698
float xBelow = -1e6, xAbove = 1e6;
666-
const auto& clus = TPCClusterIdxStruct.clusters[sector][row][icl0];
667-
float xc, yc, zc;
668-
mTPCCorrMapsLoader.Transform(sector, row, clus.getPad(), clus.getTime(), xc, yc, zc, mctr.bcInTF / 8.); // nominal time of the track
669699

670700
for (int entID = entTRefIDs.first; entID < entTRefIDs.second; entID++) {
671701
const auto& refTr = mSelTRefs[entID];
672702
if (refTr.getUserField() != sector % 18) {
673703
continue;
674704
}
675-
if (refTr.getX() < xc && refTr.getX() > xBelow && refTr.getX() > xc - params.maxTPCRefExtrap) {
705+
if ((refTr.getX() < xc) && (refTr.getX() > xBelow) && (refTr.getX() > xc - params.maxTPCRefExtrap)) {
676706
xBelow = refTr.getX();
677707
entIDBelow = entID;
678708
}
679-
if (refTr.getX() > xc && refTr.getX() < xAbove && refTr.getX() < xc + params.maxTPCRefExtrap) {
709+
if ((refTr.getX() > xc) && (refTr.getX() < xAbove) && (refTr.getX() < xc + params.maxTPCRefExtrap)) {
680710
xAbove = refTr.getX();
681711
entIDAbove = entID;
682712
}
@@ -688,42 +718,53 @@ void TrackMCStudy::fillMCClusterInfo(const o2::globaltracking::RecoContainer& re
688718
o2::track::TrackPar tparAbove, tparBelow;
689719
bool okBelow = entIDBelow >= 0 && prop->PropagateToXBxByBz((tparBelow = mSelTRefs[entIDBelow]), xc, 0.99, 2.);
690720
bool okAbove = entIDAbove >= 0 && prop->PropagateToXBxByBz((tparAbove = mSelTRefs[entIDAbove]), xc, 0.99, 2.);
691-
if ((!okBelow && !okAbove) || (params.requireTopBottomRefs && (!okBelow || !okAbove)) || (params.rejectClustersResStat > 0. && gRandom->Rndm() < params.rejectClustersResStat)) {
721+
if ((!okBelow && !okAbove) || (params.requireTopBottomRefs && (!okBelow || !okAbove))) {
692722
continue;
693723
}
694-
ClResTPC clRes{};
695724

696725
int nmeas = 0;
726+
auto& clCont = clRes.contTracks.emplace_back();
727+
clCont.corrAttach = corrAttach;
697728
if (okBelow) {
698-
clRes.below = {mSelTRefs[entIDBelow].getX(), tparBelow.getY(), tparBelow.getZ()};
699-
clRes.snp += tparBelow.getSnp();
700-
clRes.tgl += tparBelow.getTgl();
729+
clCont.below = {mSelTRefs[entIDBelow].getX(), tparBelow.getY(), tparBelow.getZ()};
730+
clCont.snp += tparBelow.getSnp();
731+
clCont.tgl += tparBelow.getTgl();
732+
clCont.q2pt += tparBelow.getQ2Pt();
701733
nmeas++;
702734
}
703735
if (okAbove) {
704-
clRes.above = {mSelTRefs[entIDAbove].getX(), tparAbove.getY(), tparAbove.getZ()};
705-
clRes.snp += tparAbove.getSnp();
706-
clRes.tgl += tparBelow.getTgl();
736+
clCont.above = {mSelTRefs[entIDAbove].getX(), tparAbove.getY(), tparAbove.getZ()};
737+
clCont.snp += tparAbove.getSnp();
738+
clCont.tgl += tparAbove.getTgl();
739+
clCont.q2pt += tparAbove.getQ2Pt();
707740
nmeas++;
708741
}
709742
if (nmeas) {
743+
if (clRes.contTracks.size() == 1) {
744+
int occBin = mctr.bcInTF / 8 * mNTPCOccBinLengthInv;
745+
clRes.occ = occBin < 0 ? mTBinClOcc[0] : (occBin >= mTBinClOcc.size() ? mTBinClOcc.back() : mTBinClOcc[occBin]);
746+
}
747+
clCont.xyz = {xc, yc, zc};
710748
if (nmeas > 1) {
711-
clRes.snp *= 0.5;
712-
clRes.tgl *= 0.5;
749+
clCont.snp *= 0.5;
750+
clCont.tgl *= 0.5;
751+
clCont.q2pt *= 0.5;
713752
}
714-
clRes.clPos = {xc, yc, zc};
715-
clRes.sect = sector;
716-
clRes.row = row;
717-
clRes.qtot = clus.getQtot();
718-
clRes.qmax = clus.getQmax();
719-
clRes.flags = clus.getFlags();
720-
clRes.ncont = ncontLb;
721-
int occBin = mctr.bcInTF / 8 * mNTPCOccBinLengthInv;
722-
clRes.occ = occBin < 0 ? mTBinClOcc[0] : (occBin >= mTBinClOcc.size() ? mTBinClOcc.back() : mTBinClOcc[occBin]);
723-
(*mDBGOut) << "clres" << "clr=" << clRes << "\n";
753+
} else {
754+
clRes.contTracks.pop_back();
724755
}
725756
}
726757
}
758+
if (clRes.getNCont()) {
759+
clRes.sect = sector;
760+
clRes.row = row;
761+
clRes.qtot = clus.getQtot();
762+
clRes.qmax = clus.getQmax();
763+
clRes.flags = clus.getFlags();
764+
clRes.ncont = ncontLb;
765+
clRes.sortCont();
766+
(*mDBGOut) << "clres" << "clr=" << clRes << "\n";
767+
}
727768
}
728769
}
729770
}

0 commit comments

Comments
 (0)