Skip to content

Commit a2c3f71

Browse files
martenolemartenole
committed
TPC SCD map creation from unbinned residuals
This allows to specify the track type(s) which should be used for the extraction of the distortion map
1 parent 085f847 commit a2c3f71

File tree

7 files changed

+151
-28
lines changed

7 files changed

+151
-28
lines changed

‎Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualReaderSpec.h‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -15,16 +15,18 @@
1515
#define O2_TPC_RESIDUALREADER
1616

1717
#include "Framework/DataProcessorSpec.h"
18+
#include "ReconstructionDataFormats/GlobalTrackID.h"
1819

1920
using namespace o2::framework;
21+
using GID = o2::dataformats::GlobalTrackID;
2022

2123
namespace o2
2224
{
2325
namespace tpc
2426
{
2527

2628
/// create a processor spec
27-
framework::DataProcessorSpec getTPCResidualReaderSpec();
29+
framework::DataProcessorSpec getTPCResidualReaderSpec(bool doBinning, GID::mask_t src);
2830

2931
} // namespace tpc
3032
} // namespace o2

‎Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -64,7 +64,8 @@ void TPCInterpolationDPL::updateTimeDependentParams(ProcessingContext& pc)
6464
bool limitTracks = (SpacePointsCalibConfParam::Instance().maxTracksPerCalibSlot < 0) ? true : false;
6565
int nTracksPerTfMax = (nTfs > 0 && !limitTracks) ? SpacePointsCalibConfParam::Instance().maxTracksPerCalibSlot / nTfs : -1;
6666
if (nTracksPerTfMax > 0) {
67-
LOGP(info, "We will stop processing tracks after validating {} tracks per TF", nTracksPerTfMax);
67+
LOGP(info, "We will stop processing tracks after validating {} tracks per TF, since we want to accumulate {} tracks for a slot with {} TFs",
68+
nTracksPerTfMax, SpacePointsCalibConfParam::Instance().maxTracksPerCalibSlot, nTfs);
6869
} else if (nTracksPerTfMax < 0) {
6970
LOG(info) << "The number of processed tracks per TF is not limited";
7071
} else {

‎Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualReaderSpec.cxx‎

Lines changed: 112 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,7 @@
2323
#include "Framework/ControlService.h"
2424
#include "DataFormatsTPC/Defs.h"
2525
#include "SpacePoints/TrackResiduals.h"
26+
#include "SpacePoints/TrackInterpolation.h"
2627
#include "DetectorsBase/Propagator.h"
2728
#include "CommonUtils/StringUtils.h"
2829
#include "TPCInterpolationWorkflow/TPCResidualReaderSpec.h"
@@ -38,7 +39,7 @@ namespace tpc
3839
class TPCResidualReader : public Task
3940
{
4041
public:
41-
TPCResidualReader() = default;
42+
TPCResidualReader(bool doBinning, GID::mask_t src) : mDoBinning(doBinning), mSources(src) {}
4243
~TPCResidualReader() override = default;
4344
void init(InitContext& ic) final;
4445
void run(ProcessingContext& pc) final;
@@ -53,10 +54,18 @@ class TPCResidualReader : public Task
5354
std::unique_ptr<TFile> mFile;
5455
std::unique_ptr<TTree> mTreeResiduals;
5556
std::unique_ptr<TTree> mTreeStats;
57+
std::unique_ptr<TTree> mTreeUnbinnedResiduals;
58+
bool mDoBinning{false};
59+
GID::mask_t mSources;
5660
TrackResiduals mTrackResiduals;
5761
std::vector<std::string> mFileNames{""}; ///< input files
5862
std::string mOutfile{"debugVoxRes.root"}; ///< output file name
59-
std::vector<o2::tpc::TrackResiduals::LocalResid> mResiduals, *mResidualsPtr = &mResiduals;
63+
std::vector<TrackResiduals::LocalResid> mResiduals, *mResidualsPtr = &mResiduals; ///< binned residuals input
64+
std::array<std::vector<TrackResiduals::LocalResid>, SECTORSPERSIDE * SIDES> mResidualsSector; ///< binned residuals generated on-the-fly
65+
std::array<std::vector<TrackResiduals::LocalResid>*, SECTORSPERSIDE * SIDES> mResidualsSectorPtr; ///< for setting branch addresses
66+
std::array<std::vector<TrackResiduals::VoxStats>, SECTORSPERSIDE * SIDES> mVoxStatsSector; ///< voxel statistics generated on-the-fly
67+
std::vector<UnbinnedResid> mUnbinnedResiduals, *mUnbinnedResidualsPtr = &mUnbinnedResiduals; ///< unbinned residuals input
68+
std::vector<TrackDataCompact> mTrackData, *mTrackDataPtr = &mTrackData; ///< the track references for unbinned residuals
6069
};
6170

6271
void TPCResidualReader::fillResiduals(const int iSec)
@@ -111,15 +120,86 @@ void TPCResidualReader::run(ProcessingContext& pc)
111120
{
112121
mTrackResiduals.createOutputFile(mOutfile.data()); // FIXME remove when map output is handled properly
113122

114-
for (int iSec = 0; iSec < SECTORSPERSIDE * SIDES; ++iSec) {
123+
if (mDoBinning) {
124+
// initialize the trees for the binned residuals and the statistics
125+
LOG(info) << "Preparing the binning of residuals";
126+
mTreeResiduals = std::make_unique<TTree>("resid", "TPC binned residuals");
127+
for (int iSec = 0; iSec < SECTORSPERSIDE * SIDES; ++iSec) {
128+
mResidualsSectorPtr[iSec] = &mResidualsSector[iSec];
129+
mVoxStatsSector[iSec].resize(mTrackResiduals.getNVoxelsPerSector());
130+
for (int ix = 0; ix < mTrackResiduals.getNXBins(); ++ix) {
131+
for (int ip = 0; ip < mTrackResiduals.getNY2XBins(); ++ip) {
132+
for (int iz = 0; iz < mTrackResiduals.getNZ2XBins(); ++iz) {
133+
auto& statsVoxel = mVoxStatsSector[iSec][mTrackResiduals.getGlbVoxBin(ix, ip, iz)];
134+
// COG estimates are set to the bin center by default
135+
mTrackResiduals.getVoxelCoordinates(iSec, ix, ip, iz, statsVoxel.meanPos[TrackResiduals::VoxX], statsVoxel.meanPos[TrackResiduals::VoxF], statsVoxel.meanPos[TrackResiduals::VoxZ]);
136+
}
137+
}
138+
}
139+
mTreeResiduals->Branch(Form("sec%d", iSec), &mResidualsSectorPtr[iSec]);
140+
}
141+
// now go through all input files, apply the track selection and fill be binned residuals
115142
for (const auto& file : mFileNames) {
116143
LOGP(info, "Processing residuals from file {}", file);
117-
118-
// set up the tree from the input file
119144
connectTree(file);
145+
for (int iEntry = 0; iEntry < mTreeUnbinnedResiduals->GetEntries(); ++iEntry) {
146+
mTreeUnbinnedResiduals->GetEntry(iEntry);
147+
for (const auto& trkInfo : mTrackData) {
148+
if (!GID::includesSource(trkInfo.sourceId, mSources)) {
149+
continue;
150+
}
151+
for (int i = trkInfo.idxFirstResidual; i < trkInfo.idxFirstResidual + trkInfo.nResiduals; ++i) {
152+
const auto& residIn = mUnbinnedResiduals[i];
153+
int sec = residIn.sec;
154+
auto& residVecOut = mResidualsSector[sec];
155+
auto& statVecOut = mVoxStatsSector[sec];
156+
std::array<unsigned char, TrackResiduals::VoxDim> bvox;
157+
float xPos = param::RowX[residIn.row];
158+
float yPos = residIn.y * param::MaxY / 0x7fff;
159+
float zPos = residIn.z * param::MaxZ / 0x7fff;
160+
if (!mTrackResiduals.findVoxelBin(sec, xPos, yPos, zPos, bvox)) {
161+
// we are not inside any voxel
162+
LOGF(debug, "Dropping residual in sec(%i), x(%f), y(%f), z(%f)", sec, xPos, yPos, zPos);
163+
continue;
164+
}
165+
residVecOut.emplace_back(residIn.dy, residIn.dz, residIn.tgSlp, bvox);
166+
auto& stat = statVecOut[mTrackResiduals.getGlbVoxBin(bvox)];
167+
float& binEntries = stat.nEntries;
168+
float oldEntries = binEntries++;
169+
float norm = 1.f / binEntries;
170+
// update COG for voxel bvox (update for X only needed in case binning is not per pad row)
171+
float xPosInv = 1.f / xPos;
172+
stat.meanPos[TrackResiduals::VoxX] = (stat.meanPos[TrackResiduals::VoxX] * oldEntries + xPos) * norm;
173+
stat.meanPos[TrackResiduals::VoxF] = (stat.meanPos[TrackResiduals::VoxF] * oldEntries + yPos * xPosInv) * norm;
174+
stat.meanPos[TrackResiduals::VoxZ] = (stat.meanPos[TrackResiduals::VoxZ] * oldEntries + zPos * xPosInv) * norm;
175+
}
176+
}
177+
}
178+
mTreeResiduals->Fill();
179+
for (auto& residuals : mResidualsSector) {
180+
residuals.clear();
181+
}
182+
}
183+
}
120184

121-
// fill the residuals for one sector
122-
fillResiduals(iSec);
185+
for (int iSec = 0; iSec < SECTORSPERSIDE * SIDES; ++iSec) {
186+
if (mDoBinning) {
187+
auto brResid = mTreeResiduals->GetBranch(Form("sec%d", iSec));
188+
brResid->SetAddress(&mResidualsPtr);
189+
190+
for (int iEntry = 0; iEntry < brResid->GetEntries(); ++iEntry) {
191+
brResid->GetEntry(iEntry);
192+
mTrackResiduals.getLocalResVec().insert(mTrackResiduals.getLocalResVec().end(), mResiduals.begin(), mResiduals.end());
193+
}
194+
mTrackResiduals.setStats(mVoxStatsSector[iSec], iSec);
195+
} else {
196+
for (const auto& file : mFileNames) {
197+
LOGP(info, "Processing residuals from file {}", file);
198+
// set up the tree from the input file
199+
connectTree(file);
200+
// fill the residuals for one sector
201+
fillResiduals(iSec);
202+
}
123203
}
124204

125205
// do processing
@@ -139,19 +219,34 @@ void TPCResidualReader::run(ProcessingContext& pc)
139219

140220
void TPCResidualReader::connectTree(const std::string& filename)
141221
{
142-
mTreeResiduals.reset(nullptr); // in case it was already loaded
143-
mTreeStats.reset(nullptr); // in case it was already loaded
222+
if (!mDoBinning) {
223+
// in case we do the binning on-the-fly, we fill these trees manually
224+
// and don't want to delete them in between
225+
mTreeResiduals.reset(nullptr); // in case it was already loaded
226+
mTreeStats.reset(nullptr); // in case it was already loaded
227+
}
228+
mTreeUnbinnedResiduals.reset(nullptr); // in case it was already loaded
144229
mFile.reset(TFile::Open(filename.c_str()));
145230
assert(mFile && !mFile->IsZombie());
146-
mTreeResiduals.reset((TTree*)mFile->Get("resid"));
147-
assert(mTreeResiduals);
148-
mTreeStats.reset((TTree*)mFile->Get("stats"));
149-
assert(mTreeStats);
150-
151-
LOG(info) << "Loaded tree from " << filename << " with " << mTreeResiduals->GetEntries() << " entries";
231+
if (!mDoBinning) {
232+
// we load the already binned residuals and the voxel statistics
233+
mTreeResiduals.reset((TTree*)mFile->Get("resid"));
234+
assert(mTreeResiduals);
235+
mTreeStats.reset((TTree*)mFile->Get("stats"));
236+
assert(mTreeStats);
237+
LOG(info) << "Loaded tree from " << filename << " with " << mTreeResiduals->GetEntries() << " entries";
238+
} else {
239+
// we load the unbinned residuals
240+
LOG(info) << "Loading the binned residuals";
241+
mTreeUnbinnedResiduals.reset((TTree*)mFile->Get("unbinnedResid"));
242+
assert(mTreeUnbinnedResiduals);
243+
mTreeUnbinnedResiduals->SetBranchAddress("res", &mUnbinnedResidualsPtr);
244+
mTreeUnbinnedResiduals->SetBranchAddress("trackInfo", &mTrackDataPtr);
245+
LOG(info) << "Loaded tree from " << filename << " with " << mTreeUnbinnedResiduals->GetEntries() << " entries";
246+
}
152247
}
153248

154-
DataProcessorSpec getTPCResidualReaderSpec()
249+
DataProcessorSpec getTPCResidualReaderSpec(bool doBinning, GID::mask_t src)
155250
{
156251
std::vector<InputSpec> inputs;
157252
std::vector<OutputSpec> outputs;
@@ -160,7 +255,7 @@ DataProcessorSpec getTPCResidualReaderSpec()
160255
"tpc-residual-reader",
161256
inputs,
162257
outputs,
163-
AlgorithmSpec{adaptFromTask<TPCResidualReader>()},
258+
AlgorithmSpec{adaptFromTask<TPCResidualReader>(doBinning, src)},
164259
Options{
165260
{"residuals-infiles", VariantType::String, "o2tpc_residuals.root", {"comma-separated list of input files or .txt file containing list of input files"}},
166261
{"input-dir", VariantType::String, "none", {"Input directory"}},

‎Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-map-creator.cxx‎

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,14 +11,18 @@
1111

1212
#include "Framework/DataProcessorSpec.h"
1313
#include "TPCInterpolationWorkflow/TPCResidualReaderSpec.h"
14+
#include "ReconstructionDataFormats/GlobalTrackID.h"
1415
#include "CommonUtils/ConfigurableParam.h"
1516

1617
using namespace o2::framework;
18+
using GID = o2::dataformats::GlobalTrackID;
1719

1820
// we need to add workflow options before including Framework/runDataProcessing
1921
void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions)
2022
{
2123
std::vector<o2::framework::ConfigParamSpec> options{
24+
{"track-sources", VariantType::String, std::string{GID::ALL}, {"comma-separated list of sources to use for tracking"}},
25+
{"start-from-unbinned", VariantType::Bool, false, {"Do the binning of the residuals on-the-fly (taking into account allowed track-sources)"}},
2226
{"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings ..."}}};
2327
std::swap(workflowOptions, options);
2428
}
@@ -31,7 +35,13 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext)
3135
{
3236
o2::conf::ConfigurableParam::updateFromString(configcontext.options().get<std::string>("configKeyValues"));
3337
o2::conf::ConfigurableParam::writeINI("o2tpcmapcreator-workflow_configuration.ini");
38+
39+
GID::mask_t allowedSources = GID::getSourcesMask("ITS-TPC,ITS-TPC-TRD,ITS-TPC-TOF,ITS-TPC-TRD-TOF");
40+
GID::mask_t src = allowedSources & GID::getSourcesMask(configcontext.options().get<std::string>("track-sources"));
41+
42+
auto doBinning = configcontext.options().get<bool>("start-from-unbinned");
43+
3444
WorkflowSpec specs;
35-
specs.emplace_back(o2::tpc::getTPCResidualReaderSpec());
45+
specs.emplace_back(o2::tpc::getTPCResidualReaderSpec(doBinning, src));
3646
return specs;
3747
}

‎Detectors/TPC/calibration/SpacePoints/include/SpacePoints/ResidualAggregator.h‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -55,16 +55,16 @@ struct ResidualsContainer {
5555
void writeToFile(bool closeFileAfterwards);
5656

5757
const TrackResiduals* trackResiduals{nullptr};
58-
std::array<std::vector<TrackResiduals::LocalResid>, SECTORSPERSIDE * SIDES> residuals{}; ///< local residuals per sector which are sent to the aggregator
58+
std::array<std::vector<TrackResiduals::LocalResid>, SECTORSPERSIDE * SIDES> residuals{}; ///< local (binned) residuals per sector
5959
std::array<std::vector<TrackResiduals::LocalResid>*, SECTORSPERSIDE * SIDES> residualsPtr{};
60-
std::array<std::vector<TrackResiduals::VoxStats>, SECTORSPERSIDE * SIDES> stats{}; ///< voxel statistics sent to the aggregator
60+
std::array<std::vector<TrackResiduals::VoxStats>, SECTORSPERSIDE * SIDES> stats{}; ///< voxel statistics per sector
6161
std::array<std::vector<TrackResiduals::VoxStats>*, SECTORSPERSIDE * SIDES> statsPtr{};
6262
std::vector<uint32_t> tfOrbits, *tfOrbitsPtr{&tfOrbits}; ///< first TF orbit
6363
std::vector<uint32_t> sumOfResiduals, *sumOfResidualsPtr{&sumOfResiduals}; ///< sum of residuals for each TF
6464
std::vector<o2::ctp::LumiInfo> lumi, *lumiPtr{&lumi}; ///< luminosity information from CTP per TF
65-
std::vector<UnbinnedResid> unbinnedRes, *unbinnedResPtr{&unbinnedRes}; ///< unbinned residuals
65+
std::vector<UnbinnedResid> unbinnedRes, *unbinnedResPtr{&unbinnedRes}; ///< unbinned residuals which are sent to the aggregator
6666
std::vector<TrackData> trkData, *trkDataPtr{&trkData}; ///< track data and cluster ranges
67-
std::vector<TrackDataCompact> trackInfo, *trackInfoPtr{&trackInfo}; ///< allows to obtain track type for each binned residual downstream
67+
std::vector<TrackDataCompact> trackInfo, *trackInfoPtr{&trackInfo}; ///< allows to obtain track type for each unbinned residual downstream
6868

6969
std::string fileName{"o2tpc_residuals"};
7070
std::unique_ptr<TFile> fileOut{nullptr};

‎Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackResiduals.h‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -420,14 +420,17 @@ class TrackResiduals
420420
/// Closes the file with the debug output.
421421
void closeOutputFile();
422422

423+
/// Set the voxel statistics directly from outside
424+
void setStats(const std::vector<TrackResiduals::VoxStats>& statsIn, int iSec);
425+
423426
/// Fill statistics from TTree
424427
void fillStats(int iSec);
425428

426429
/// clear member to be able to process new sector or new input files
427430
void clear();
428431

429432
private:
430-
bool mInitResultsContainer{false};
433+
std::bitset<SECTORSPERSIDE * SIDES> mInitResultsContainer{};
431434

432435
// some constants
433436
static constexpr float sFloatEps{1.e-7f}; ///< float epsilon for robust linear fitting

‎Detectors/TPC/calibration/SpacePoints/src/TrackResiduals.cxx‎

Lines changed: 15 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -173,10 +173,10 @@ void TrackResiduals::initBinning()
173173
//______________________________________________________________________________
174174
void TrackResiduals::initResultsContainer(int iSec)
175175
{
176-
if (mInitResultsContainer) {
176+
if (mInitResultsContainer.test(iSec)) {
177177
return;
178178
}
179-
mInitResultsContainer = true;
179+
mInitResultsContainer.set(iSec);
180180
mVoxelResults[iSec].resize(mNVoxPerSector);
181181
for (int ix = 0; ix < mNXBins; ++ix) {
182182
for (int ip = 0; ip < mNY2XBins; ++ip) {
@@ -300,6 +300,18 @@ void TrackResiduals::setKernelType(KernelType kernel, float bwX, float bwP, floa
300300
///
301301
///////////////////////////////////////////////////////////////////////////////
302302

303+
void TrackResiduals::setStats(const std::vector<TrackResiduals::VoxStats>& statsIn, int iSec)
304+
{
305+
initResultsContainer(iSec);
306+
std::vector<VoxRes>& secDataTmp = mVoxelResults[iSec];
307+
for (int iVox = 0; iVox < mNVoxPerSector; ++iVox) {
308+
secDataTmp[iVox].stat[VoxX] = statsIn[iVox].meanPos[VoxX];
309+
secDataTmp[iVox].stat[VoxF] = statsIn[iVox].meanPos[VoxF];
310+
secDataTmp[iVox].stat[VoxZ] = statsIn[iVox].meanPos[VoxZ];
311+
secDataTmp[iVox].stat[VoxDim] = statsIn[iVox].nEntries;
312+
}
313+
}
314+
303315
void TrackResiduals::fillStats(int iSec)
304316
{
305317
initResultsContainer(iSec);
@@ -1452,5 +1464,5 @@ void TrackResiduals::printMem() const
14521464
void TrackResiduals::clear()
14531465
{
14541466
getLocalResVec().clear();
1455-
mInitResultsContainer = false;
1467+
mInitResultsContainer.reset();
14561468
}

0 commit comments

Comments
 (0)