Skip to content

Commit 6b9b5b9

Browse files
martenolemartenole
committed
Macro for creating average distortion maps
1 parent 0b5d4db commit 6b9b5b9

2 files changed

Lines changed: 315 additions & 0 deletions

File tree

‎Detectors/TPC/calibration/SpacePoints/CMakeLists.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -34,3 +34,7 @@ o2_target_root_dictionary(SpacePoints
3434
include/SpacePoints/ResidualAggregator.h
3535
include/SpacePoints/SpacePointsCalibConfParam.h
3636
LINKDEF src/SpacePointCalibLinkDef.h)
37+
38+
o2_add_test_root_macro(macro/staticMapCreator.C
39+
PUBLIC_LINK_LIBRARIES O2::SpacePoints
40+
LABELS tpc COMPILE_ONLY)

‎Detectors/TPC/calibration/SpacePoints/macro/staticMapCreator.C‎

Lines changed: 311 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,311 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
#if !defined(__CLING__) || defined(__ROOTCLING__)
13+
#include "CCDB/CcdbApi.h"
14+
#include "CCDB/BasicCCDBManager.h"
15+
16+
#include "Framework/Logger.h"
17+
#include "CommonConstants/LHCConstants.h"
18+
#include "SpacePoints/SpacePointsCalibConfParam.h"
19+
#include "SpacePoints/TrackResiduals.h"
20+
#include "SpacePoints/TrackInterpolation.h"
21+
#include "DataFormatsParameters/GRPMagField.h"
22+
#include "DataFormatsTPC/Defs.h"
23+
#include "ReconstructionDataFormats/GlobalTrackID.h"
24+
#include "DetectorsBase/MatLayerCylSet.h"
25+
#include "DetectorsBase/Propagator.h"
26+
27+
#include <TFile.h>
28+
#include <TTree.h>
29+
#include <TGeoManager.h>
30+
#include <TGrid.h>
31+
32+
#include <fstream>
33+
#include <string>
34+
#include <vector>
35+
#include <memory>
36+
#include <array>
37+
38+
#include <boost/algorithm/string/predicate.hpp>
39+
#include <boost/filesystem.hpp>
40+
41+
#else
42+
43+
#error This macro must run in compiled mode
44+
45+
#endif
46+
47+
using namespace o2::tpc;
48+
using GID = o2::dataformats::GlobalTrackID;
49+
50+
std::vector<string> getInputFileList(const std::string& fileInput)
51+
{
52+
std::vector<std::string> fileList;
53+
std::vector<std::string> fileListVerified;
54+
// check if only one input file (a txt file contaning a list of files is provided)
55+
if (boost::algorithm::ends_with(fileInput, "txt")) {
56+
LOGP(info, "Reading files from input file list {}", fileInput);
57+
std::ifstream is(fileInput);
58+
std::istream_iterator<std::string> start(is);
59+
std::istream_iterator<std::string> end;
60+
fileList.insert(fileList.begin(), start, end);
61+
} else {
62+
fileList.push_back(fileInput);
63+
}
64+
65+
for (const auto& file : fileList) {
66+
if ((file.find("alien://") == 0) && !gGrid && !TGrid::Connect("alien://")) {
67+
LOG(fatal) << "Failed to open alien connection";
68+
}
69+
std::unique_ptr<TFile> filePtr(TFile::Open(file.data()));
70+
if (!filePtr || !filePtr->IsOpen() || filePtr->IsZombie()) {
71+
LOGP(warning, "Could not open file {}", file);
72+
continue;
73+
}
74+
fileListVerified.push_back(file);
75+
}
76+
77+
if (fileListVerified.size() == 0) {
78+
LOGP(error, "No input files to process");
79+
}
80+
return fileListVerified;
81+
}
82+
83+
bool revalidateTrack(const TrackData& trk, const SpacePointsCalibConfParam& params)
84+
{
85+
if (trk.nClsITS < params.minITSNCls) {
86+
return false;
87+
}
88+
if (trk.nClsTPC < params.minTPCNCls) {
89+
return false;
90+
}
91+
if (trk.nTrkltsTRD > 0 && trk.nTrkltsTRD < params.minTRDNTrklts) {
92+
// in case nTrkltsTRD == 0 this is an ITS-TPC-TOF track which we might not want to cut
93+
return false;
94+
}
95+
// track quality cuts
96+
if (trk.chi2ITS / trk.nClsITS > params.maxITSChi2) {
97+
return false;
98+
}
99+
if (trk.chi2TPC / trk.nClsTPC > params.maxTPCChi2) {
100+
return false;
101+
}
102+
if (trk.nTrkltsTRD > 0 && trk.chi2TRD / trk.nTrkltsTRD > params.maxTRDChi2) {
103+
return false;
104+
}
105+
106+
if (params.cutOnDCA) {
107+
auto propagator = o2::base::Propagator::Instance();
108+
o2::track::TrackPar trkPar(trk.x, trk.alpha, trk.p);
109+
// o2::track::TrackPar trkPar = trk.par; // for the next version of o2::tpc::TrackData where its stored as TrackPar directly
110+
if (!propagator->propagateToX(trkPar, 0, propagator->getNominalBz())) {
111+
return false;
112+
}
113+
if (trkPar.getX() * trkPar.getX() + trkPar.getY() * trkPar.getY() > params.maxDCA * params.maxDCA) {
114+
LOGP(debug, "DCA cut not passed {}", std::sqrt(trkPar.getX() * trkPar.getX() + trkPar.getY() * trkPar.getY()));
115+
return false;
116+
}
117+
LOGP(debug, "DCA cut OK {}", std::sqrt(trkPar.getX() * trkPar.getX() + trkPar.getY() * trkPar.getY()));
118+
}
119+
return true;
120+
}
121+
122+
void staticMapCreator(std::string fileInput = "files.txt",
123+
int runNumber = 527976,
124+
std::string fileOutput = "voxRes.root",
125+
std::string trackSources = static_cast<std::string>(GID::ALL))
126+
{
127+
128+
// Obtain configuration
129+
const SpacePointsCalibConfParam& params = SpacePointsCalibConfParam::Instance();
130+
if (!boost::filesystem::exists("scdconfig.ini")) {
131+
LOG(warn) << "Did not find configuration file. Using default parameters and storing them in scdconfig.ini";
132+
params.writeINI("scdconfig.ini", "scdcalib"); // to write default parameters to a file
133+
} else {
134+
params.updateFromFile("scdconfig.ini");
135+
}
136+
LOG(info) << "----- Dumping configuration values START -----";
137+
params.printKeyValues();
138+
LOG(info) << "----- Dumping configuration values END -----";
139+
140+
GID::mask_t allowedSources = GID::getSourcesMask("ITS-TPC,ITS-TPC-TRD,ITS-TPC-TOF,ITS-TPC-TRD-TOF");
141+
GID::mask_t sources = allowedSources & GID::getSourcesMask(trackSources);
142+
143+
// Get CCDB objects
144+
auto& ccdbmgr = o2::ccdb::BasicCCDBManager::instance();
145+
ccdbmgr.setURL("https://alice-ccdb.cern.ch");
146+
auto runDuration = ccdbmgr.getRunDuration(runNumber);
147+
auto tRun = runDuration.first + (runDuration.second - runDuration.first) / 2; // time stamp for the middle of the run duration
148+
ccdbmgr.setTimestamp(tRun);
149+
150+
// CTP orbit reset time
151+
auto orbitResetTimeNS = ccdbmgr.get<std::vector<int64_t>>("CTP/Calib/OrbitReset");
152+
int64_t orbitResetTimeMS = (*orbitResetTimeNS)[0] * 1e-3;
153+
LOGP(info, "Orbit reset time in MS is {}", orbitResetTimeMS);
154+
155+
auto geoAligned = ccdbmgr.get<TGeoManager>("GLO/Config/GeometryAligned");
156+
auto magField = ccdbmgr.get<o2::parameters::GRPMagField>("GLO/Config/GRPMagField");
157+
const o2::base::MatLayerCylSet* matLut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdbmgr.get<o2::base::MatLayerCylSet>("GLO/Param/MatLUT"));
158+
o2::base::Propagator::initFieldFromGRP(magField);
159+
auto prop = o2::base::Propagator::Instance();
160+
prop->setMatLUT(matLut);
161+
162+
// Input
163+
auto fileList = getInputFileList(fileInput);
164+
165+
std::array<std::vector<TrackResiduals::LocalResid>, SECTORSPERSIDE * SIDES> binnedResidualsSec; // binned residuals generated on-the-fly
166+
std::array<std::vector<TrackResiduals::LocalResid>*, SECTORSPERSIDE * SIDES> binnedResidualsSecPtr; // for setting branch addresses
167+
std::array<std::vector<TrackResiduals::VoxStats>, SECTORSPERSIDE * SIDES> voxStatsSec; // voxel statistics generated on-the-fly
168+
std::vector<TrackResiduals::LocalResid> binnedResiduals, *binnedResidualsPtr = &binnedResiduals; // binned residuals
169+
170+
TrackResiduals trackResiduals;
171+
trackResiduals.init();
172+
trackResiduals.createOutputFile(fileOutput.c_str());
173+
174+
std::unique_ptr<TTree> treeBinnedResiduals = std::make_unique<TTree>("resid", "TPC binned residuals");
175+
if (!params.writeBinnedResiduals) {
176+
treeBinnedResiduals->SetDirectory(nullptr);
177+
}
178+
for (int iSec = 0; iSec < SECTORSPERSIDE * SIDES; ++iSec) {
179+
binnedResidualsSecPtr[iSec] = &binnedResidualsSec[iSec];
180+
voxStatsSec[iSec].resize(trackResiduals.getNVoxelsPerSector());
181+
for (int ix = 0; ix < trackResiduals.getNXBins(); ++ix) {
182+
for (int ip = 0; ip < trackResiduals.getNY2XBins(); ++ip) {
183+
for (int iz = 0; iz < trackResiduals.getNZ2XBins(); ++iz) {
184+
auto& statsVoxel = voxStatsSec[iSec][trackResiduals.getGlbVoxBin(ix, ip, iz)];
185+
// COG estimates are set to the bin center by default
186+
trackResiduals.getVoxelCoordinates(iSec, ix, ip, iz, statsVoxel.meanPos[TrackResiduals::VoxX], statsVoxel.meanPos[TrackResiduals::VoxF], statsVoxel.meanPos[TrackResiduals::VoxZ]);
187+
}
188+
}
189+
}
190+
treeBinnedResiduals->Branch(Form("sec%d", iSec), &binnedResidualsSecPtr[iSec]);
191+
}
192+
193+
std::unique_ptr<TFile> inputFile;
194+
std::unique_ptr<TTree> treeUnbinnedResiduals;
195+
std::unique_ptr<TTree> treeTrackData;
196+
std::unique_ptr<TTree> treeRecords;
197+
std::vector<UnbinnedResid> unbinnedResiduals, *unbinnedResidualsPtr = &unbinnedResiduals; // unbinned residuals input
198+
std::vector<TrackDataCompact> trackRefs, *trackRefsPtr = &trackRefs; // the track references for unbinned residuals
199+
std::vector<TrackData> trackData, *trackDataPtr = &trackData; // additional track information (chi2, nClusters, track parameters)
200+
std::vector<uint32_t> orbits, *orbitsPtr = &orbits; // first orbit for each TF in the input data
201+
202+
for (const auto& fileName : fileList) {
203+
LOGP(info, "Processing input file {}", fileName);
204+
treeUnbinnedResiduals.reset(nullptr);
205+
treeTrackData.reset(nullptr);
206+
treeRecords.reset(nullptr);
207+
inputFile.reset(TFile::Open(fileName.c_str()));
208+
if (!inputFile || inputFile->IsZombie()) {
209+
LOGP(info, "Skipping file {}", fileName);
210+
continue;
211+
}
212+
treeUnbinnedResiduals.reset((TTree*)inputFile->Get("unbinnedResid"));
213+
treeUnbinnedResiduals->SetBranchAddress("res", &unbinnedResidualsPtr);
214+
treeUnbinnedResiduals->SetBranchAddress("trackInfo", &trackRefsPtr);
215+
if (params.useTrackData) {
216+
treeTrackData.reset((TTree*)inputFile->Get("trackData"));
217+
treeTrackData->SetBranchAddress("trk", &trackDataPtr);
218+
if (treeTrackData->GetEntries() != treeUnbinnedResiduals->GetEntries()) {
219+
LOGP(error, "The input trees with unbinned residuals and track information have a different number of entries ({} vs {})",
220+
treeUnbinnedResiduals->GetEntries(), treeTrackData->GetEntries());
221+
}
222+
}
223+
treeRecords.reset((TTree*)inputFile->Get("records"));
224+
treeRecords->SetBranchAddress("firstTForbit", &orbitsPtr);
225+
treeRecords->GetEntry(0); // per input file there is only a single entry in the tree
226+
for (int iEntry = 0; iEntry < treeUnbinnedResiduals->GetEntries(); ++iEntry) {
227+
if (params.timeFilter) {
228+
int64_t tfTimeInMS = orbitResetTimeMS + orbits[iEntry] * o2::constants::lhc::LHCOrbitMUS * 1.e-3;
229+
if (tfTimeInMS < params.startTimeMS || tfTimeInMS > params.endTimeMS) {
230+
LOGP(debug, "Dropping TF at index {} with time {} and orbit {}", iEntry, tfTimeInMS, orbits[iEntry]);
231+
continue;
232+
}
233+
}
234+
treeUnbinnedResiduals->GetEntry(iEntry);
235+
if (params.useTrackData) {
236+
treeTrackData->GetEntry(iEntry);
237+
}
238+
auto nTracks = trackRefs.size();
239+
for (size_t iTrack = 0; iTrack < nTracks; ++iTrack) {
240+
const auto& trkInfo = trackRefs[iTrack];
241+
if (!GID::includesSource(trkInfo.sourceId, sources)) {
242+
continue;
243+
}
244+
if (params.useTrackData) {
245+
const auto& trk = trackData[iTrack];
246+
if (!revalidateTrack(trk, params)) {
247+
continue;
248+
}
249+
}
250+
for (unsigned int i = trkInfo.idxFirstResidual; i < trkInfo.idxFirstResidual + trkInfo.nResiduals; ++i) {
251+
const auto& residIn = unbinnedResiduals[i];
252+
int sec = residIn.sec;
253+
auto& residVecOut = binnedResidualsSec[sec];
254+
auto& statVecOut = voxStatsSec[sec];
255+
std::array<unsigned char, TrackResiduals::VoxDim> bvox;
256+
float xPos = param::RowX[residIn.row];
257+
float yPos = residIn.y * param::MaxY / 0x7fff;
258+
float zPos = residIn.z * param::MaxZ / 0x7fff;
259+
if (!trackResiduals.findVoxelBin(sec, xPos, yPos, zPos, bvox)) {
260+
// we are not inside any voxel
261+
LOGF(debug, "Dropping residual in sec(%i), x(%f), y(%f), z(%f)", sec, xPos, yPos, zPos);
262+
continue;
263+
}
264+
residVecOut.emplace_back(residIn.dy, residIn.dz, residIn.tgSlp, bvox);
265+
auto& stat = statVecOut[trackResiduals.getGlbVoxBin(bvox)];
266+
float& binEntries = stat.nEntries;
267+
float oldEntries = binEntries++;
268+
float norm = 1.f / binEntries;
269+
// update COG for voxel bvox (update for X only needed in case binning is not per pad row)
270+
float xPosInv = 1.f / xPos;
271+
stat.meanPos[TrackResiduals::VoxX] = (stat.meanPos[TrackResiduals::VoxX] * oldEntries + xPos) * norm;
272+
stat.meanPos[TrackResiduals::VoxF] = (stat.meanPos[TrackResiduals::VoxF] * oldEntries + yPos * xPosInv) * norm;
273+
stat.meanPos[TrackResiduals::VoxZ] = (stat.meanPos[TrackResiduals::VoxZ] * oldEntries + zPos * xPosInv) * norm;
274+
}
275+
}
276+
}
277+
treeBinnedResiduals->Fill();
278+
for (auto& resid : binnedResidualsSec) {
279+
resid.clear();
280+
}
281+
}
282+
283+
for (int iSec = 0; iSec < SECTORSPERSIDE * SIDES; ++iSec) {
284+
// for each sector fill the vector of local residuals from the respective branch
285+
auto brResid = treeBinnedResiduals->GetBranch(Form("sec%d", iSec));
286+
brResid->SetAddress(&binnedResidualsPtr);
287+
for (int iEntry = 0; iEntry < brResid->GetEntries(); ++iEntry) {
288+
brResid->GetEntry(iEntry);
289+
trackResiduals.getLocalResVec().insert(trackResiduals.getLocalResVec().end(), binnedResiduals.begin(), binnedResiduals.end());
290+
}
291+
trackResiduals.setStats(voxStatsSec[iSec], iSec);
292+
// do processing
293+
trackResiduals.processSectorResiduals(iSec);
294+
// do cleanup
295+
trackResiduals.clear();
296+
}
297+
298+
if (params.writeBinnedResiduals) {
299+
trackResiduals.getOutputFilePtr()->cd();
300+
treeBinnedResiduals->Write();
301+
}
302+
treeBinnedResiduals.reset();
303+
trackResiduals.closeOutputFile();
304+
305+
treeUnbinnedResiduals.reset();
306+
treeTrackData.reset();
307+
treeRecords.reset();
308+
inputFile.reset();
309+
310+
LOG(info) << "Done processing";
311+
}

0 commit comments

Comments
 (0)