From b2a33732afa4685ea78b0c9075912f2e6c3a3a31 Mon Sep 17 00:00:00 2001 From: Ole Schmidt Date: Mon, 20 Apr 2020 10:52:05 +0200 Subject: [PATCH 1/4] TOF DPL reader for matching information --- .../tofworkflow/CMakeLists.txt | 5 +- .../TOFWorkflow/TOFMatchedReaderSpec.h | 58 +++++++++++ .../tofworkflow/src/TOFMatchedReaderSpec.cxx | 99 +++++++++++++++++++ 3 files changed, 160 insertions(+), 2 deletions(-) create mode 100644 Detectors/GlobalTrackingWorkflow/tofworkflow/include/TOFWorkflow/TOFMatchedReaderSpec.h create mode 100644 Detectors/GlobalTrackingWorkflow/tofworkflow/src/TOFMatchedReaderSpec.cxx diff --git a/Detectors/GlobalTrackingWorkflow/tofworkflow/CMakeLists.txt b/Detectors/GlobalTrackingWorkflow/tofworkflow/CMakeLists.txt index 631a8603c4594..0e2aed78c0327 100644 --- a/Detectors/GlobalTrackingWorkflow/tofworkflow/CMakeLists.txt +++ b/Detectors/GlobalTrackingWorkflow/tofworkflow/CMakeLists.txt @@ -10,8 +10,9 @@ o2_add_library(TOFWorkflow SOURCES src/RecoWorkflowSpec.cxx - src/TOFMatchedWriterSpec.cxx - src/TOFCalibWriterSpec.cxx + src/TOFMatchedWriterSpec.cxx + src/TOFCalibWriterSpec.cxx + src/TOFMatchedReaderSpec.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::TOFBase O2::DataFormatsTOF O2::TOFReconstruction O2::GlobalTracking O2::GlobalTrackingWorkflow O2::TOFWorkflowUtils O2::TOFCalibration O2::DataFormatsFT0 O2::FT0Reconstruction O2::FITWorkflow) diff --git a/Detectors/GlobalTrackingWorkflow/tofworkflow/include/TOFWorkflow/TOFMatchedReaderSpec.h b/Detectors/GlobalTrackingWorkflow/tofworkflow/include/TOFWorkflow/TOFMatchedReaderSpec.h new file mode 100644 index 0000000000000..d8b1488a57d9f --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/tofworkflow/include/TOFWorkflow/TOFMatchedReaderSpec.h @@ -0,0 +1,58 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @file TOFMatchedReaderSpec.h + +#ifndef O2_TOF_MATCHINFOREADER +#define O2_TOF_MATCHINFOREADER + +#include "TFile.h" +#include "TTree.h" + +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" +#include "ReconstructionDataFormats/MatchInfoTOF.h" +#include "SimulationDataFormat/MCCompLabel.h" + +namespace o2 +{ +namespace tof +{ + +class TOFMatchedReader : public o2::framework::Task +{ + public: + TOFMatchedReader(bool useMC) : mUseMC(useMC) {} + ~TOFMatchedReader() override = default; + void init(o2::framework::InitContext& ic) final; + void run(o2::framework::ProcessingContext& pc) final; + + private: + void connectTree(const std::string& filename); + + bool mUseMC = false; + std::string mInFileName{"o2match_tof.root"}; + std::string mInTreeName{"matchTOF"}; + std::unique_ptr mFile = nullptr; + std::unique_ptr mTree = nullptr; + std::vector mMatches, *mMatchesPtr = &mMatches; + std::vector mLabelTOF, *mLabelTOFPtr = &mLabelTOF; + std::vector mLabelTPC, *mLabelTPCPtr = &mLabelTPC; + std::vector mLabelITS, *mLabelITSPtr = &mLabelITS; +}; + +/// create a processor spec +/// read matched TOF clusters from a ROOT file +framework::DataProcessorSpec getTOFMatchedReaderSpec(bool useMC); + +} // namespace tof +} // namespace o2 + +#endif /* O2_TOF_MATCHINFOREADER */ diff --git a/Detectors/GlobalTrackingWorkflow/tofworkflow/src/TOFMatchedReaderSpec.cxx b/Detectors/GlobalTrackingWorkflow/tofworkflow/src/TOFMatchedReaderSpec.cxx new file mode 100644 index 0000000000000..93c432cfa21ae --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/tofworkflow/src/TOFMatchedReaderSpec.cxx @@ -0,0 +1,99 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @file TOFMatchedReaderSpec.cxx + +#include + +#include "TTree.h" +#include "TFile.h" + +#include "TOFWorkflow/TOFMatchedReaderSpec.h" +#include "Framework/ControlService.h" +#include "Framework/ConfigParamRegistry.h" +#include "Headers/DataHeader.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "ReconstructionDataFormats/MatchInfoTOF.h" + +using namespace o2::framework; + +namespace o2 +{ +namespace tof +{ + +void TOFMatchedReader::init(InitContext& ic) +{ + // get the option from the init context + LOG(INFO) << "Init TOF matching info reader!"; + mInFileName = ic.options().get("tof-matched-infile"); + mInTreeName = ic.options().get("treename"); + connectTree(mInFileName); +} + +void TOFMatchedReader::connectTree(const std::string& filename) +{ + mTree.reset(nullptr); // in case it was already loaded + mFile.reset(TFile::Open(filename.c_str())); + assert(mFile && !mFile->IsZombie()); + mTree.reset((TTree*)mFile->Get(mInTreeName.c_str())); + assert(mTree); + mTree->SetBranchAddress("TOFMatchInfo", &mMatchesPtr); + if (mUseMC) { + mTree->SetBranchAddress("MatchTOFMCTruth", &mLabelTOFPtr); + mTree->SetBranchAddress("MatchTPCMCTruth", &mLabelTPCPtr); + mTree->SetBranchAddress("MatchITSMCTruth", &mLabelITSPtr); + } + LOG(INFO) << "Loaded tree from " << filename << " with " << mTree->GetEntries() << " entries"; +} + +void TOFMatchedReader::run(ProcessingContext& pc) +{ + auto currEntry = mTree->GetReadEntry() + 1; + assert(currEntry < mTree->GetEntries()); // this should not happen + mTree->GetEntry(currEntry); + LOG(INFO) << "Pushing " << mMatches.size() << " TOF matchings at entry " << currEntry; + + pc.outputs().snapshot(Output{o2::header::gDataOriginTOF, "MATCHINFOS", 0, Lifetime::Timeframe}, mMatches); + if (mUseMC) { + pc.outputs().snapshot(Output{o2::header::gDataOriginTOF, "MATCHTOFINFOSMC", 0, Lifetime::Timeframe}, mLabelTOF); + pc.outputs().snapshot(Output{o2::header::gDataOriginTOF, "MATCHTPCINFOSMC", 0, Lifetime::Timeframe}, mLabelTPC); + pc.outputs().snapshot(Output{o2::header::gDataOriginTOF, "MATCHITSINFOSMC", 0, Lifetime::Timeframe}, mLabelITS); + } + + if (mTree->GetReadEntry() + 1 >= mTree->GetEntries()) { + pc.services().get().endOfStream(); + pc.services().get().readyToQuit(QuitRequest::Me); + } +} + +DataProcessorSpec getTOFMatchedReaderSpec(bool useMC) +{ + std::vector outputs; + outputs.emplace_back(o2::header::gDataOriginTOF, "MATCHINFOS", 0, Lifetime::Timeframe); + if (useMC) { + outputs.emplace_back(o2::header::gDataOriginTOF, "MATCHTOFINFOSMC", 0, Lifetime::Timeframe); + outputs.emplace_back(o2::header::gDataOriginTOF, "MATCHTPCINFOSMC", 0, Lifetime::Timeframe); + outputs.emplace_back(o2::header::gDataOriginTOF, "MATCHITSINFOSMC", 0, Lifetime::Timeframe); + } + + return DataProcessorSpec{ + "TOFMatchedReader", + Inputs{}, + outputs, + AlgorithmSpec{adaptFromTask(useMC)}, + Options{ + {"tof-matched-infile", VariantType::String, "o2match_tof.root", {"Name of the input file"}}, + {"treename", VariantType::String, "matchTOF", {"Name of top-level TTree"}}, + }}; +} +} // namespace tof +} // namespace o2 From 4d9374fee01118f5e1cbabb22f0ed296e88a38cb Mon Sep 17 00:00:00 2001 From: Ole Schmidt Date: Mon, 20 Apr 2020 10:52:55 +0200 Subject: [PATCH 2/4] Fix short range moving average filter and add possibility to dump outlier information for thesis plots --- .../SpacePoints/SpacePointsCalibParam.h | 2 +- .../include/SpacePoints/TrackResiduals.h | 32 +++++++--- .../SpacePoints/src/TrackResiduals.cxx | 61 +++++++++++++++---- 3 files changed, 75 insertions(+), 20 deletions(-) diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibParam.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibParam.h index 28151063be032..48ad4eb80cb4b 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibParam.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibParam.h @@ -16,7 +16,7 @@ #ifndef ALICEO2_TPC_PARAM_H_ #define ALICEO2_TPC_PARAM_H_ -#define TPC_RUN2 +//#define TPC_RUN2 #include "DataFormatsTPC/Constants.h" diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackResiduals.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackResiduals.h index 3ad2a98c10ab0..5be72b9045e5e 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackResiduals.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackResiduals.h @@ -115,6 +115,17 @@ class TrackResiduals UShort_t npValid{0}; ///< number of valid TPC clusters }; + struct DebugOutliers { + int idx{-1}; + std::array x{}; + std::array xFlagged{}; + std::array dY{}; + std::array dZ{}; + std::array residHelixY{}; + std::array residHelixZ{}; + unsigned char flags{0}; + }; + // -------------------------------------- initialization -------------------------------------------------- /// Steers the initialization (binning, default settings for smoothing, container for the results). void init(); @@ -275,7 +286,7 @@ class TrackResiduals /// \param res[1] contains the offset (b) bool fitPoly1(int nCl, std::array& x, std::array& y, std::array& res); - /// For a given set of points, calculate the differences from each point to the fitted lines from all other points in their neighbourhoods (+- mNMALong points) + /// For a given set of points, calculate the differences from each point to the fitted lines from all other points in their neighbourhoods (+- mNMAShort points) void diffToLocLine(int np, int idxOffset, const std::array& x, const std::array& y, std::array& diffY); /// For a given set of points, calculate their deviation from the moving average (build from the neighbourhood +- mNMALong points) @@ -424,6 +435,7 @@ class TrackResiduals // -------------------------------------- settings -------------------------------------------------- + void setPathToResFileRun2(std::string fPath) { mPathToResidualFiles = fPath; } void setLocalResFileName(std::string fName) { mLocalResFileName = fName; } void setLocalResTreeName(std::string tName) { mLocalResTreeName = tName; } void setLocalResBranchName(std::string bName) { mLocalResBranchName = bName; } @@ -477,9 +489,14 @@ class TrackResiduals /// \param fName Filename void dumpToFile(const std::vector& vec, const std::string fName) const; - /// Dumps the content of a vector to output tree (mTreeOut) to compare it to the data in the AliRoot version - /// \param vec Data vector with floating point values - void dumpVector(const std::vector& vec); + /// Dumps the content of an array to the specified file (used for visualization of outlier rejection) + /// \param arr Data array + /// \param fName Filename + void dumpArrayToFile(const std::array& arr, const std::string fName) const; + + /// Dumps the collected data from tracks with information from the outlier rejection routines + /// \param vec Data vector with all relevant information on a per track basis + void dumpTracks(const std::vector& vec); /// Dumps the full results for a given sector to the debug tree (only if an output file has been created before). /// \param iSec Sector to dump @@ -548,7 +565,7 @@ class TrackResiduals int mMinEntriesPerVoxel{15}; ///< minimum number of points in voxel for processing float mLTMCut{.75f}; ///< fraction op points to keep when trimming input data float mMinFracLTM{.5f}; ///< minimum fraction of points to keep when trimming data to fit expected sigma - float mMinValidVoxFracDrift{.5f}; ///< if more than this fraction of bins are bad for one pad row the bad row is declared bad + float mMinValidVoxFracDrift{.5f}; ///< if more than this fraction of bins are bad for one pad row the whole pad row is declared bad int mMinGoodXBinsToCover{3}; ///< minimum number of consecutive good bins, otherwise bins are declared bad int mMaxBadXBinsToCover{4}; ///< a lower number of consecutive bad X bins will not be declared bad float mMaxFracBadRowsPerSector{.4f}; ///< maximum fraction of bad rows before whole sector is masked @@ -581,6 +598,7 @@ class TrackResiduals std::unique_ptr mRun2DeltaTree{}; ///< tree with Run 2 cluster residuals bool mFilterOutliers = true; ///< flag, if outliers from the cluster residual trees should be rejected int mNMALong{15}; ///< number of points to be used for moving average (long range) + int mNMAShort{3}; ///< number of points to be used for estimation of distance from local line (short range) float mMaxRejFrac{.15f}; ///< if the fraction of rejected clusters of a track is higher, the full track is invalidated float mMaxRMSLong{.8f}; ///< maximum variance of the cluster residuals wrt moving avarage for a track to be considered // buffer arrays as in AliTPCDcalibRes @@ -599,8 +617,8 @@ class TrackResiduals float mTgl{0.f}; ///< fitted track dip angle int mNCl{0}; ///< number of clusters in the track // debugging - std::vector mOutVector{}; ///< this vector can be filled with data to stream it to a ROOT tree - std::vector* mOutVectorPtr{&mOutVector}; ///< pointer to set the branch address of the debug ROOT tree to + std::vector mOutVector{}; ///< this vector can be filled with data to stream it to a ROOT tree + std::vector* mOutVectorPtr{&mOutVector}; ///< pointer to set the branch address of the debug ROOT tree to }; //_____________________________________________________ diff --git a/Detectors/TPC/calibration/SpacePoints/src/TrackResiduals.cxx b/Detectors/TPC/calibration/SpacePoints/src/TrackResiduals.cxx index d7787169c30e4..dc610f2a9efe7 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/TrackResiduals.cxx +++ b/Detectors/TPC/calibration/SpacePoints/src/TrackResiduals.cxx @@ -38,7 +38,7 @@ #include -#define TPC_RUN2 // if defined, use run 2 geometry for TPC +//#define TPC_RUN2 // if defined, use run 2 geometry for TPC #define LOCAL_RESIDUAL_FORMAT_OLD // if defined, data in compact trees is stored as Double32_t, otherwise as short #ifdef LOCAL_RESIDUAL_FORMAT_OLD @@ -356,6 +356,7 @@ void TrackResiduals::buildLocalResidualTreesFromRun2Data() std::array counterTrkValidation{0}; int nRejCl = 0, nRejHelix = 0, nRejQpt = 0, nRejValidation = 0; LOG(info) << "Building local residual trees from " << nTracks << " tracks."; + std::vector debugOutliers; for (int iTrk = 0; iTrk < nTracks; ++iTrk) { /* printf("Checking track %i\n", iTrk); @@ -363,6 +364,7 @@ void TrackResiduals::buildLocalResidualTreesFromRun2Data() break; } */ + DebugOutliers debug; brTRDOK->GetEntry(iTrk); brITSOK->GetEntry(iTrk); if (!mDeltaStruct.trdOK || !mDeltaStruct.itsOK) { @@ -413,20 +415,33 @@ void TrackResiduals::buildLocalResidualTreesFromRun2Data() ++nTracksSelectedWithOutliers; //printf("Checking track %i\n", iTrk); bool resHelix = compareToHelix(residHelixY, residHelixZ); + debug.idx = iTrk; + debug.x = mArrR; + debug.dY = mArrDY; + debug.dZ = mArrDZ; + debug.residHelixY = residHelixY; + debug.residHelixZ = residHelixZ; + /* printf("Printing helix residuals for track %i\n", iTrk); for (int i = 0; i < param::NPadRows; ++i) { printf("residHelixY[%03i]=% .4f \t \t residHelixZ[%03i]=% .4f\n", i, residHelixY[i], i, residHelixZ[i]); } */ + if (mFilterOutliers && !resHelix) { // too strong deviation to helix -> discard track + //printf("Track %i deviates strongly from helix\n", iTrk); ++nRejHelix; + debug.flags = 1 << 0; // rejected by helix fit + debugOutliers.push_back(debug); continue; } if (fabsf(mQpt) > param::MaxQ2Pt) { // discard low pt tracks now that a more precise q/pt estimate is available ++nRejQpt; + debug.flags = 1 << 1; // rejected by qpt cut + debugOutliers.push_back(debug); continue; } @@ -479,10 +494,18 @@ void TrackResiduals::buildLocalResidualTreesFromRun2Data() // done converting everything to sector frame ++mNCl; } + debug.x = mArrX; + debug.dY = mArrDY; + debug.dZ = mArrDZ; if (mFilterOutliers && !validateTrack(counterTrkValidation)) { ++nRejValidation; + debug.flags = 1 << 2; + debugOutliers.push_back(debug); continue; } + debug.xFlagged = mArrX; + debug.flags = 1 << 3; // good track + debugOutliers.push_back(debug); ++nTracksSelected; fillLocalResidualsTrees(); @@ -490,6 +513,7 @@ void TrackResiduals::buildLocalResidualTreesFromRun2Data() printf("Rejected due to Nclusters(%i), HelixFit(%i), qpt(%i), validation(%i)\n", nRejCl, nRejHelix, nRejQpt, nRejValidation); printf("validation failed %i times because of fraction of rej. cls and %i times because of rms and %i rest\n", counterTrkValidation[1], counterTrkValidation[2], counterTrkValidation[0]); LOG(info) << "Accepted " << nTracksSelected << " tracks. With outliers it would be " << nTracksSelectedWithOutliers; + dumpTracks(debugOutliers); writeLocalResidualTreesToFile(); } @@ -622,7 +646,6 @@ int TrackResiduals::checkResiduals(std::bitset& rejCl, float& r float average = 0.f; float rms = 0.f; for (int i = 0; i < nAcc; ++i) { - // what about points without enough neighbours?? don't they distort the average? average += yDiffLong[i]; rms += yDiffLong[i] * yDiffLong[i]; } @@ -920,8 +943,8 @@ void TrackResiduals::processSectorResiduals(int iSec) } else { smooth(iSec); } - dumpResults(iSec); - return; + //dumpResults(iSec); + //return; // process dispersions dyVec.clear(); @@ -1609,6 +1632,7 @@ void TrackResiduals::diffToMA(int np, const std::array& } int nPoints = iRight - iLeft; if (nPoints < mNMALong) { + // this cannot happen, since at least mNMALong points are required as neighbours for this function to be called continue; } float movingAverage = (sum[iRight] - sum[iLeft - 1] - (sum[i] - sum[i - 1])) / nPoints; @@ -1646,8 +1670,8 @@ void TrackResiduals::diffToLocLine(int np, int idxOffset, const std::array& residHel void TrackResiduals::fitCircle(int nCl, std::array& x, std::array& y, float& xc, float& yc, float& r, std::array& residHelixY) { + // this fast algebraic circle fit is described here: + // https://dtcenter.org/met/users/docs/write_ups/circle_fit.pdf float xMean = 0.f; float yMean = 0.f; @@ -2145,6 +2171,17 @@ bool TrackResiduals::fitPoly1(int nCl, std::array& x, st /// /////////////////////////////////////////////////////////////////////////////// +void TrackResiduals::dumpArrayToFile(const std::array& arr, const std::string fName) const +{ + std::ofstream fOut(fName.data()); + if (fOut.is_open()) { + for (const auto& elem : arr) { + fOut << std::fixed << std::setprecision(6) << elem << std::endl; + } + fOut.close(); + } +} + void TrackResiduals::dumpToFile(const std::vector& vec, const std::string fName = "output.txt") const { std::ofstream fOut(fName.data()); @@ -2158,10 +2195,10 @@ void TrackResiduals::dumpToFile(const std::vector& vec, const std::string void TrackResiduals::createOutputFile() { - mFileOut = std::make_unique("voxelResultsO2.root", "recreate"); - mTreeOut = std::make_unique("debugTree", "voxel results"); - mTreeOut->Branch("voxRes", &mVoxelResultsOutPtr); - //mTreeOut->Branch("debug", &mOutVectorPtr); + mFileOut = std::make_unique("debugOutliers.root", "recreate"); + mTreeOut = std::make_unique("debugTree", "outliers"); + //mTreeOut->Branch("voxRes", &mVoxelResultsOutPtr); + mTreeOut->Branch("debug", &mOutVectorPtr); } void TrackResiduals::closeOutputFile() @@ -2173,7 +2210,7 @@ void TrackResiduals::closeOutputFile() mFileOut.reset(); } -void TrackResiduals::dumpVector(const std::vector& vec) +void TrackResiduals::dumpTracks(const std::vector& vec) { if (mTreeOut) { mOutVector = vec; From bc841034cdde3660d5d6d373cda10142b9f208a3 Mon Sep 17 00:00:00 2001 From: Ole Schmidt Date: Mon, 20 Apr 2020 10:53:56 +0200 Subject: [PATCH 3/4] Prepare TPC track interpolation for DPL driven input and suppress tree input for TPC SCD TrackInterpolation --- .../include/SpacePoints/TrackInterpolation.h | 119 ++++------- .../SpacePoints/src/SpacePointCalibLinkDef.h | 2 + .../SpacePoints/src/TrackInterpolation.cxx | 193 +++++------------- 3 files changed, 90 insertions(+), 224 deletions(-) diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h index 66934e7ed4978..1cb4c61e29880 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h @@ -16,6 +16,7 @@ #ifndef ALICEO2_TPC_TRACKINTERPOLATION_H_ #define ALICEO2_TPC_TRACKINTERPOLATION_H_ +#include #include "CommonDataFormat/EvIndex.h" #include "CommonDataFormat/RangeReference.h" #include "ReconstructionDataFormats/Track.h" @@ -23,7 +24,6 @@ #include "ReconstructionDataFormats/MatchInfoTOF.h" #include "DataFormatsITSMFT/Cluster.h" #include "DataFormatsITS/TrackITS.h" -#include "DataFormatsTPC/ClusterNativeHelper.h" #include "DataFormatsTPC/TrackTPC.h" #include "DataFormatsTPC/Constants.h" #include "DataFormatsTOF/Cluster.h" @@ -54,6 +54,7 @@ struct TPCClusterResiduals { void setZ(float val) { z = fabs(val) < param::MaxZ ? static_cast(val * 0x7fff / param::MaxZ) : static_cast(std::copysign(1., val) * 0x7fff); } void setPhi(float val) { phi = fabs(val) < param::MaxTgSlp ? static_cast(val * 0x7fff / param::MaxTgSlp) : static_cast(std::copysign(1., val) * 0x7fff); } void setTgl(float val) { tgl = fabs(val) < param::MaxTgSlp ? static_cast(val * 0x7fff / param::MaxTgSlp) : static_cast(std::copysign(1., val) * 0x7fff); } + ClassDefNV(TPCClusterResiduals, 1); }; /// Structure filled for each track with track quality information and a vector with TPCClusterResiduals @@ -70,6 +71,7 @@ struct TrackData { // TODO: Add an additional structure with event information and reference to the tracks for given event? int nTracksInEvent{}; ///< total number of tracks in given event / timeframe ? Used for multiplicity estimate o2::dataformats::RangeReference<> clIdx{}; ///< index of first cluster residual and total number of cluster residuals of this track + ClassDefNV(TrackData, 1); }; /// \class TrackInterpolation @@ -105,6 +107,8 @@ class TrackInterpolation unsigned short clAvailable{0}; }; + // -------------------------------------- processing functions -------------------------------------------------- + /// Initialize everything void init(); @@ -122,7 +126,7 @@ class TrackInterpolation /// Interpolate ITS-TPC-TOF tracks in the TPC and store residuals to TPC clusters /// \param matchTOF /// \return flag if track could be processed successfully - bool interpolateTrackITSTOF(const std::pair, o2::dataformats::MatchInfoTOF>& matchTOF); + bool interpolateTrackITSTOF(const o2::dataformats::MatchInfoTOF& matchTOF); /// Check if given ITS-TPC track fullfills quality criteria /// \param matchITSTPC ITS-TPC track to be checked @@ -130,41 +134,12 @@ class TrackInterpolation /// \return Flag if quality criteria are met bool trackPassesQualityCuts(const o2::dataformats::TrackTPCITS& matchITSTPC, bool hasOuterPoint = true) const; - /// Load specified entry for all input trees - /// \param iEntry number of tree entry - void loadEntryForTrees(int iEntry); - - /// Attach all input trees - void attachInputTrees(); - - /// Add branches to the output trees - void prepareOutputTrees(); /// Reset cache and output vectors void reset(); // -------------------------------------- settings -------------------------------------------------- - /// Sets the tree/chain containing ITS-TPC matched tracks - void setInputTreeITSTPCTracks(TTree* tree) { mTreeITSTPCTracks = tree; } - /// Sets the tree/chain containing TPC tracks - void setInputTreeTPCTracks(TTree* tree) { mTreeTPCTracks = tree; } - ///< Sets the reader for TPC clusters - void setInputTPCClustersReader(ClusterNativeHelper::Reader* reader) { mTPCClusterReader = reader; } - /// Sets the tree/chain containing ITS tracks - void setInputTreeITSTracks(TTree* tree) { mTreeITSTracks = tree; } - /// Sets the tree/chain containing ITS clusters - void setInputTreeITSClusters(TTree* tree) { mTreeITSClusters = tree; } - /// Sets the tree/chain containing TOF matching information - void setInputTreeTOFMatches(TTree* tree) { mTreeTOFMatches = tree; } - ///< Sets the tree/chain containing TOF clusters - void setInputTreeTOFClusters(TTree* tree) { mTreeTOFClusters = tree; } - - /// Sets the output tree for track information - void setOutputTreeTracks(TTree* tree) { mTreeOutTrackData = tree; } - /// Sets the output tree for TPC cluster residuals - void setOutputTreeResiduals(TTree* tree) { mTreeOutClusterRes = tree; } - /// Sets the maximum phi angle at which track extrapolation is aborted void setMaxSnp(float snp) { mMaxSnp = snp; } /// Sets the maximum step length for track extrapolation @@ -174,80 +149,56 @@ class TrackInterpolation /// Sets whether ITS tracks without match in TRD or TOF should be processed as well void setDoITSOnlyTracks(bool flag) { mDoITSOnlyTracks = flag; } - /// Setters and Getters for the input branch/file names - void setITSTPCTrackBranchName(const std::string& name) { mITSTPCTrackBranchName = name; } - void setTPCTrackBranchName(const std::string& name) { mTPCTrackBranchName = name; } - void setTPCTrackClIdxBranchName(const std::string& name) { mTPCTrackClIdxBranchName = name; } - void setTPCClusterFileName(const std::string& name) { mTPCClusterFileName = name; } - void setITSTrackBranchName(const std::string& name) { mITSTrackBranchName = name; } - void setITSClusterBranchName(const std::string& name) { mITSClusterBranchName = name; } - void setTOFMatchingBranchName(const std::string& name) { mTOFMatchingBranchName = name; } - void setTOFClusterBranchName(const std::string& name) { mTOFClusterBranchName = name; } - - const std::string& getITSTPCTrackBranchName() const { return mITSTPCTrackBranchName; } - const std::string& getTPCTrackBranchName() const { return mTPCTrackBranchName; } - const std::string& getTPCTrackClIdxBranchName() const { return mTPCTrackClIdxBranchName; } - const std::string& getTPCClusterFileName() const { return mTPCClusterFileName; } - const std::string& getITSTrackBranchName() const { return mITSTrackBranchName; } - const std::string& getITSClusterBranchName() const { return mITSClusterBranchName; } - const std::string& getTOFMatchingBranchName() const { return mTOFMatchingBranchName; } - const std::string& getTOFClusterBranchName() const { return mTOFClusterBranchName; } + // --------------------------------- input --------------------------------------------- + + /// Sets input ITS tracks received via DPL + void setITSTracksInp(const gsl::span inp) { mITSTracksArray = inp; } + /// Sets input ITS-TPC matched tracks received via DPL + void setITSTPCTrackMatchesInp(const gsl::span inp) { mITSTPCTracksArray = inp; } + /// Sets input TPC tracks received via DPL + void setTPCTracksInp(const gsl::span inp) { mTPCTracksArray = inp; } + /// Sets input TPC track cluster indices received via DPL + void setTPCTrackClusIdxInp(const gsl::span inp) { mTPCTracksClusIdx = inp; } + /// Sets input TPC clusters received via DPL + void setTPCClustersInp(const o2::tpc::ClusterNativeAccess* inp) { mTPCClusterIdxStruct = inp; } + /// Sets input TOF matching records received via DPL + void setTOFMatchesInp(const gsl::span inp) { mTOFMatchesArray = inp; } + /// Sets input TOF clusters received via DPL + void setTOFClustersInp(const gsl::span inp) { mTOFClustersArray = inp; } + + // --------------------------------- output --------------------------------------------- + std::vector& getClusterResiduals() { return mClRes; } + std::vector& getReferenceTracks() { return mTrackData; } private: - // names of input files / branches - std::string mITSTPCTrackBranchName{"TPCITS"}; ///< name of branch containing ITS-TPC matched tracks - std::string mTPCTrackBranchName{"Tracks"}; ///< name of branch containing TPC tracks (needed for TPC cluster access) - std::string mTPCTrackClIdxBranchName{"ClusRefs"}; ///< name of branch containing TPC tracks cluster refs (needed for TPC cluster access) - std::string mTPCClusterFileName{"tpc-native-clusters.root"}; ///< name of file containing TPC native clusters - std::string mITSTrackBranchName{"ITSTrack"}; ///< name of branch containing input ITS tracks - std::string mITSClusterBranchName{"ITSCluster"}; ///< name of branch containing input ITS clusters - std::string mTOFMatchingBranchName{"TOFMatchInfo"}; ///< name of branch containing matching info for ITS-TPC-TOF tracks - std::string mTOFClusterBranchName{"TOFCluster"}; ///< name of branch containing TOF clusters - - // parameters + // parameters + settings float mTPCTimeBinMUS{.2f}; ///< TPC time bin duration in us float mSigYZ2TOF{.75f}; ///< for now assume cluster error for TOF equal for all clusters in both Y and Z - - // settings float mMaxSnp{.85f}; ///< max snp when propagating ITS tracks float mMaxStep{2.f}; ///< maximum step for propagation - int mMatCorr{1}; ///< if material correction should be done + int mMatCorr{0}; ///< if material correction should be done bool mDoITSOnlyTracks{false}; ///< if ITS only tracks should be processed or not // input - TTree* mTreeITSTPCTracks{nullptr}; ///< input tree for ITS-TPC matched tracks - TTree* mTreeTPCTracks{nullptr}; ///< input tree for TPC tracks - TTree* mTreeITSTracks{nullptr}; ///< input tree for ITS tracks - TTree* mTreeITSClusters{nullptr}; ///< input tree for ITS clusters - TTree* mTreeTOFMatches{nullptr}; ///< input tree for ITS-TPC-TOF matches - TTree* mTreeTOFClusters{nullptr}; ///< input tree for TOF clusters - std::vector* mITSTPCTrackVecInput{nullptr}; ///< input vector with ITS-TPC matched tracks - std::vector* mTPCTrackVecInput{nullptr}; ///< input vector with TPC tracks - std::vector* mTPCTrackClIdxVecInput{nullptr}; ///< input vector with TPC tracks cluster indicies - std::vector* mITSTrackVecInput{nullptr}; ///< input vector with ITS tracks - std::vector* mITSClusterVecInput{nullptr}; ///< input vector with ITS clusters - ClusterNativeHelper::Reader* mTPCClusterReader = nullptr; ///< TPC cluster reader - std::unique_ptr mTPCClusterIdxStructOwn; ///< struct holding the TPC cluster indices (for tree-based I/O) - std::unique_ptr mTPCClusterBufferOwn; ///< buffer for clusters in mTPCClusterIdxStructOwn - MCLabelContainer mTPCClusterMCBufferOwn; ///< buffer for mc labels - std::vector, o2::dataformats::MatchInfoTOF>>* mTOFMatchVecInput{nullptr}; ///< input vector for ITS-TPC-TOF matching information - std::vector* mTOFClusterVecInput{nullptr}; ///< input vector with TOF clusters + gsl::span mITSTPCTracksArray; ///< input ITS-TPC matched tracks from span + gsl::span mTPCTracksArray; ///< input TPC tracks from span + gsl::span mTPCTracksClusIdx; ///< input TPC tracks cluster indicies from span + gsl::span mITSTracksArray; ///< input ITS tracks from span + gsl::span mTOFMatchesArray; ///< input ITS-TPC-TOF matching information from span + gsl::span mTOFClustersArray; ///< input TOF clusters from span const ClusterNativeAccess* mTPCClusterIdxStruct = nullptr; ///< struct holding the TPC cluster indices // output - TTree* mTreeOutTrackData{nullptr}; ///< output tree with track quality information std::vector mTrackData{}; ///< this vector is used to store the track quality information on a per track basis - std::vector* mTrackDataPtr{&mTrackData}; ///< pointer to vector with track data - TTree* mTreeOutClusterRes{nullptr}; ///< output tree with TPC cluster residuals std::vector mClRes{}; ///< residuals for each available TPC cluster of all tracks - std::vector* mClResPtr{&mClRes}; ///< pointer to vector with residuals // cache std::array mCache{{}}; ///< caching positions, covariances and angles for track extrapolations and interpolation // helpers std::unique_ptr mFastTransform{}; ///< TPC cluster transformation + bool mInitDone{false}; ///< initialization done flag }; } // namespace tpc diff --git a/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h b/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h index 7961f3bc770ef..2ed783553e4ef 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h +++ b/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h @@ -17,7 +17,9 @@ #pragma link C++ class o2::tpc::TrackInterpolation + ; #pragma link C++ class o2::tpc::TrackResiduals + ; #pragma link C++ class o2::tpc::TrackData + ; +#pragma link C++ class std::vector < o2::tpc::TrackData> + ; #pragma link C++ class o2::tpc::TPCClusterResiduals + ; +#pragma link C++ class std::vector < o2::tpc::TPCClusterResiduals> + ; #pragma link C++ class o2::tpc::TrackResiduals::LocalResid + ; #pragma link C++ class o2::tpc::TrackResiduals::VoxRes + ; diff --git a/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx b/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx index b0dacba1036a5..22fb6a57d53e9 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx +++ b/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx @@ -26,7 +26,11 @@ using namespace o2::tpc; void TrackInterpolation::init() { // perform initialization - attachInputTrees(); + LOG(INFO) << "Start initializing TrackInterpolation"; + if (mInitDone) { + LOG(error) << "Initialization already performed."; + return; + } const auto& elParam = ParameterElectronics::Instance(); mTPCTimeBinMUS = elParam.ZbinWidth; @@ -34,60 +38,70 @@ void TrackInterpolation::init() std::unique_ptr fastTransform = (TPCFastTransformHelperO2::instance()->create(0)); mFastTransform = std::move(fastTransform); - if (mTreeOutTrackData && mTreeOutClusterRes) { - prepareOutputTrees(); - } + mInitDone = true; + LOG(INFO) << "Done initializing TrackInterpolation"; } void TrackInterpolation::process() { // main processing function - loadEntryForTrees(0); // TODO for now all input data is stored in trees with a single entry + if (!mInitDone) { + LOG(error) << "Initialization not yet done. Aborting..."; + return; + } + +#ifdef TPC_RUN2 + // processing will not work if the run 2 geometry is defined in the parameter class SpacePointsCalibParam.h + LOG(FATAL) << "Run 2 parameters compiled for the TPC geometry. Creating residual trees from Run 3 data will not work. Aborting..."; + return; +#endif std::set tracksDone; // to store indices of ITS-TPC matched tracks that have been processed - LOG(info) << "Processing " << mTOFMatchVecInput->size() << " ITS-TPC-TOF matched tracks out of " - << mITSTPCTrackVecInput->size() << " ITS-TPC matched tracks"; + LOG(INFO) << "Processing " << mTOFMatchesArray.size() << " ITS-TPC-TOF matched tracks out of " + << mITSTPCTracksArray.size() << " ITS-TPC matched tracks"; // TODO reserve only a fraction of the needed space for all tracks? How many tracks pass on average the quality cuts with how many TPC clusters? - mTrackData.reserve(mTOFMatchVecInput->size()); - mClRes.reserve(mTOFMatchVecInput->size() * param::NPadRows); + mTrackData.reserve(mTOFMatchesArray.size()); + mClRes.reserve(mTOFMatchesArray.size() * param::NPadRows); - int nTracksTPC = mTPCTrackVecInput->size(); + int nTracksTPC = mTPCTracksArray.size(); - for (const auto& trkTOF : *mTOFMatchVecInput) { + for (const auto& trkTOF : mTOFMatchesArray) { // process ITS-TPC-TOF matched tracks - if (!trackPassesQualityCuts(mITSTPCTrackVecInput->at(trkTOF.first.getIndex()))) { + if (!trackPassesQualityCuts(mITSTPCTracksArray[trkTOF.getTrackIndex()])) { + LOG(DEBUG) << "Abandoning track due to bad quality"; continue; } mTrackData.emplace_back(); if (!interpolateTrackITSTOF(trkTOF)) { + LOG(DEBUG) << "Failed to interpolate ITS-TOF track"; mTrackData.pop_back(); continue; } mTrackData.back().nTracksInEvent = nTracksTPC; - tracksDone.insert(trkTOF.first.getIndex()); + tracksDone.insert(trkTOF.getTrackIndex()); } - LOG(info) << "Could process " << tracksDone.size() << " ITS-TPC-TOF matched tracks successfully"; + LOG(INFO) << "Could process " << tracksDone.size() << " ITS-TPC-TOF matched tracks successfully"; if (mDoITSOnlyTracks) { size_t nTracksDoneITS = 0; size_t nTracksSkipped = 0; - for (std::size_t iTrk = 0; iTrk < mITSTPCTrackVecInput->size(); ++iTrk) { + for (std::size_t iTrk = 0; iTrk < mITSTPCTracksArray.size(); ++iTrk) { // process ITS-TPC matched tracks that were not matched to TOF if (tracksDone.find(iTrk) != tracksDone.end()) { // track also has a matching cluster in TOF and has already been processed ++nTracksSkipped; continue; } - const auto& trk = mITSTPCTrackVecInput->at(iTrk); + const auto& trk = mITSTPCTracksArray[iTrk]; if (!trackPassesQualityCuts(trk, false)) { continue; } mTrackData.emplace_back(); - const auto& trkTPC = mTPCTrackVecInput->at(trk.getRefTPC()); - const auto& trkITS = mITSTrackVecInput->at(trk.getRefITS()); + const auto& trkTPC = mTPCTracksArray[trk.getRefTPC()]; + const auto& trkITS = mITSTracksArray[trk.getRefITS()]; if (!extrapolateTrackITS(trkITS, trkTPC, trk.getTimeMUS().getTimeStamp(), trk.getRefTPC())) { mTrackData.pop_back(); continue; @@ -95,31 +109,26 @@ void TrackInterpolation::process() mTrackData.back().nTracksInEvent = nTracksTPC; ++nTracksDoneITS; } - LOG(info) << "Could process " << nTracksDoneITS << " ITS-TPC matched tracks successfully"; - LOG(info) << "Skipped " << nTracksSkipped << " tracks, as they were successfully propagated to TOF"; - } - - if (mTreeOutTrackData) { - mTreeOutTrackData->Fill(); - } - if (mTreeOutClusterRes) { - mTreeOutClusterRes->Fill(); + LOG(INFO) << "Could process " << nTracksDoneITS << " ITS-TPC matched tracks successfully"; + LOG(INFO) << "Skipped " << nTracksSkipped << " tracks, as they were successfully propagated to TOF"; } } bool TrackInterpolation::trackPassesQualityCuts(const o2::dataformats::TrackTPCITS& matchITSTPC, bool hasOuterPoint) const { // apply track quality cuts (assume different settings for track with and without points in TRD or TOF) - const auto& trkTPC = mTPCTrackVecInput->at(matchITSTPC.getRefTPC()); - const auto& trkITS = mITSTrackVecInput->at(matchITSTPC.getRefITS()); + const auto& trkTPC = mTPCTracksArray[matchITSTPC.getRefTPC()]; + const auto& trkITS = mITSTracksArray[matchITSTPC.getRefITS()]; if (hasOuterPoint) { // track has a match in TRD or TOF if (trkTPC.getNClusterReferences() < param::MinTPCNCls || trkITS.getNumberOfClusters() < param::MinITSNCls) { + printf("TPC clusters (%i), ITS clusters(%i)\n", trkTPC.getNClusterReferences(), trkITS.getNumberOfClusters()); return false; } if (trkTPC.getChi2() / trkTPC.getNClusterReferences() > param::MaxTPCChi2 || trkITS.getChi2() / trkITS.getNumberOfClusters() > param::MaxITSChi2) { + printf("TPC reduced chi2 (%.2f), ITS reduced chi2 (%.2f)\n", trkTPC.getChi2() / trkTPC.getNClusterReferences(), trkITS.getChi2() / trkITS.getNumberOfClusters()); return false; } } else { @@ -136,18 +145,18 @@ bool TrackInterpolation::trackPassesQualityCuts(const o2::dataformats::TrackTPCI return true; } -bool TrackInterpolation::interpolateTrackITSTOF(const std::pair, o2::dataformats::MatchInfoTOF>& matchTOF) +bool TrackInterpolation::interpolateTrackITSTOF(const o2::dataformats::MatchInfoTOF& matchTOF) { // get TPC cluster residuals to ITS-TOF only tracks size_t trkIdx = mTrackData.size() - 1; auto propagator = o2::base::Propagator::Instance(); - const auto& matchITSTPC = mITSTPCTrackVecInput->at(matchTOF.first.getIndex()); - const auto& clTOF = mTOFClusterVecInput->at(matchTOF.second.getTOFClIndex()); + const auto& matchITSTPC = mITSTPCTracksArray[matchTOF.getTrackIndex()]; + const auto& clTOF = mTOFClustersArray[matchTOF.getTOFClIndex()]; //const int clTOFSec = (TMath::ATan2(-clTOF.getY(), -clTOF.getX()) + o2::constants::math::PI) * o2::constants::math::Rad2Deg * 0.05; // taken from TOF cluster class as there is no const getter for the sector const int clTOFSec = clTOF.getCount(); const float clTOFAlpha = o2::utils::Sector2Angle(clTOFSec); - const auto& trkTPC = mTPCTrackVecInput->at(matchITSTPC.getRefTPC()); - const auto& trkITS = mITSTrackVecInput->at(matchITSTPC.getRefITS()); + const auto& trkTPC = mTPCTracksArray[matchITSTPC.getRefTPC()]; + const auto& trkITS = mITSTracksArray[matchITSTPC.getRefITS()]; auto trkWork = trkITS.getParamOut(); // reset the cache array (sufficient to set ) for (auto& elem : mCache) { @@ -161,8 +170,7 @@ bool TrackInterpolation::interpolateTrackITSTOF(const std::pair clTPCYZ; mFastTransform->TransformIdeal(sector, row, clTPC.getPad(), clTPC.getTime(), clTPCX, clTPCYZ[0], clTPCYZ[1], clusterTimeBinOffset); @@ -179,11 +187,11 @@ bool TrackInterpolation::interpolateTrackITSTOF(const std::pairPropagateToXBxByBz(trkWork, param::RowX[iRow], o2::constants::physics::MassPionCharged, mMaxSnp, mMaxStep, mMatCorr)) { - LOG(debug) << "Failed on first extrapolation"; + LOG(DEBUG) << "Failed on first extrapolation"; return false; } mCache[iRow].y[ExtOut] = trkWork.getY(); @@ -198,7 +206,7 @@ bool TrackInterpolation::interpolateTrackITSTOF(const std::pair clTOFYZ{clTOF.getY(), clTOF.getZ()}; std::array clTOFCov{mSigYZ2TOF, 0.f, mSigYZ2TOF}; // assume no correlation between y and z and equal cluster error sigma^2 = (3cm)^2 / 12 if (!propagator->PropagateToXBxByBz(trkWork, clTOFX, o2::constants::physics::MassPionCharged, mMaxSnp, mMaxStep, mMatCorr)) { - LOG(debug) << "Failed final propagation to TOF radius"; + LOG(DEBUG) << "Failed final propagation to TOF radius"; return false; } if (!trkWork.update(clTOFYZ, clTOFCov)) { - LOG(debug) << "Failed to update extrapolated ITS track with TOF cluster"; + LOG(DEBUG) << "Failed to update extrapolated ITS track with TOF cluster"; return false; } @@ -223,11 +231,11 @@ bool TrackInterpolation::interpolateTrackITSTOF(const std::pairPropagateToXBxByBz(trkWork, param::RowX[iRow], o2::constants::physics::MassPionCharged, mMaxSnp, mMaxStep, mMatCorr)) { - LOG(debug) << "Failed on back propagation"; + LOG(DEBUG) << "Failed on back propagation"; //printf("trkX(%.2f), clX(%.2f), clY(%.2f), clZ(%.2f), alphaTOF(%.2f)\n", trkWork.getX(), param::RowX[iRow], clTOFYZ[0], clTOFYZ[1], clTOFAlpha); return false; } @@ -282,6 +290,7 @@ bool TrackInterpolation::interpolateTrackITSTOF(const std::pairTransformIdeal(sector, row, cl.getPad(), cl.getTime(), x, y, z, clusterTimeBinOffset); if (!trk.rotate(o2::utils::Sector2Angle(sector))) { @@ -335,101 +343,6 @@ bool TrackInterpolation::extrapolateTrackITS(const o2::its::TrackITS& trkITS, co return true; } -void TrackInterpolation::loadEntryForTrees(int iEntry) -{ - mTreeITSTPCTracks->GetEntry(iEntry); - mTreeTPCTracks->GetEntry(iEntry); - mTreeITSTracks->GetEntry(iEntry); - mTreeITSClusters->GetEntry(iEntry); - mTreeTOFMatches->GetEntry(iEntry); - mTreeTOFClusters->GetEntry(iEntry); - mTPCClusterReader->read(iEntry); - mTPCClusterReader->fillIndex(*mTPCClusterIdxStructOwn.get(), mTPCClusterBufferOwn, mTPCClusterMCBufferOwn); - mTPCClusterIdxStruct = mTPCClusterIdxStructOwn.get(); -} - -void TrackInterpolation::attachInputTrees() -{ - // access input for ITS-TPC matched tracks (needed for the references to ITS/TPC tracks) - if (!mTreeITSTPCTracks) { - LOG(fatal) << "The input tree for ITS-TPC matched tracks is not set!"; - } - if (!mTreeITSTPCTracks->GetBranch(mITSTPCTrackBranchName.data())) { - LOG(fatal) << "Did not find ITS-TPC matched tracks branch " << mITSTPCTrackBranchName << " in the input tree"; - } - mTreeITSTPCTracks->SetBranchAddress(mITSTPCTrackBranchName.data(), &mITSTPCTrackVecInput); - LOG(info) << "Attached ITS-TPC tracks " << mITSTPCTrackBranchName << " branch with " - << mTreeITSTPCTracks->GetEntries() << " entries"; - // access input for TPC tracks (only used for cluster access) - if (!mTreeTPCTracks) { - LOG(fatal) << "The input tree for TPC tracks is not set!"; - } - if (!mTreeTPCTracks->GetBranch(mTPCTrackBranchName.data())) { - LOG(fatal) << "Did not find TPC tracks branch " << mTPCTrackBranchName << " in the input tree"; - } - mTreeTPCTracks->SetBranchAddress(mTPCTrackBranchName.data(), &mTPCTrackVecInput); - LOG(info) << "Attached TPC tracks " << mTPCTrackBranchName << " branch with " - << mTreeTPCTracks->GetEntries() << " entries"; - mTreeTPCTracks->SetBranchAddress(mTPCTrackClIdxBranchName.data(), &mTPCTrackClIdxVecInput); - LOG(info) << "Attached TPC tracks cluster indices " << mTPCTrackClIdxBranchName << " branch with " - << mTreeTPCTracks->GetBranch(mTPCTrackClIdxBranchName.data())->GetEntries() << " entries"; - - // access input for TPC clusters - if (!mTPCClusterReader) { - LOG(fatal) << "TPC clusters reader is not set"; - } - LOG(info) << "Attached TPC clusters reader with " << mTPCClusterReader->getTreeSize() << "entries"; - mTPCClusterIdxStructOwn = std::make_unique(); - // access input for ITS tracks - if (!mTreeITSTracks) { - LOG(fatal) << "The input tree for ITS tracks is not set!"; - } - if (!mTreeITSTracks->GetBranch(mITSTrackBranchName.data())) { - LOG(fatal) << "Did not find ITS tracks branch " << mITSTrackBranchName << " in the input tree"; - } - mTreeITSTracks->SetBranchAddress(mITSTrackBranchName.data(), &mITSTrackVecInput); - LOG(info) << "Attached ITS tracks " << mITSTrackBranchName << " branch with " - << mTreeITSTracks->GetEntries() << " entries"; - // access input for ITS clusters - if (!mTreeITSClusters) { - LOG(fatal) << "The input tree for ITS clusters is not set!"; - } - if (!mTreeITSClusters->GetBranch(mITSClusterBranchName.data())) { - LOG(fatal) << "Did not find ITS clusters branch " << mITSClusterBranchName << " in the input tree"; - } - mTreeITSClusters->SetBranchAddress(mITSClusterBranchName.data(), &mITSClusterVecInput); - LOG(info) << "Attached ITS clusters " << mITSClusterBranchName << " branch with " - << mTreeITSClusters->GetEntries() << " entries"; - // access input for TPC-TOF matching information - if (!mTreeTOFMatches) { - LOG(fatal) << "The input tree for with TOF matching information is not set!"; - } - if (!mTreeTOFMatches->GetBranch(mTOFMatchingBranchName.data())) { - LOG(fatal) << "Did not find TOF matches branch " << mTOFMatchingBranchName << " in the input tree"; - } - mTreeTOFMatches->SetBranchAddress(mTOFMatchingBranchName.data(), &mTOFMatchVecInput); - LOG(info) << "Attached TOF matches " << mTOFMatchingBranchName << " branch with " - << mTreeTOFMatches->GetEntries() << " entries"; - // access input for TOF clusters - if (!mTreeTOFClusters) { - LOG(fatal) << "The input tree for with TOF clusters is not set!"; - } - if (!mTreeTOFClusters->GetBranch(mTOFClusterBranchName.data())) { - LOG(fatal) << "Did not find TOF clusters branch " << mTOFClusterBranchName << " in the input tree"; - } - mTreeTOFClusters->SetBranchAddress(mTOFClusterBranchName.data(), &mTOFClusterVecInput); - LOG(info) << "Attached TOF clusters " << mTOFClusterBranchName << " branch with " - << mTreeTOFClusters->GetEntries() << " entries"; - - // TODO: add MC information -} - -void TrackInterpolation::prepareOutputTrees() -{ - mTreeOutTrackData->Branch("tracks", &mTrackDataPtr); - mTreeOutClusterRes->Branch("residuals", &mClResPtr); -} - void TrackInterpolation::reset() { mTrackData.clear(); From 464bd975bccf87b886b7c57b79b16525d46ca4f8 Mon Sep 17 00:00:00 2001 From: Ole Schmidt Date: Mon, 20 Apr 2020 10:54:46 +0200 Subject: [PATCH 4/4] Add possibility to run TPC track interpolation as DPL workflow --- .../GlobalTrackingWorkflow/CMakeLists.txt | 2 +- .../tpcinterpolationworkflow/CMakeLists.txt | 28 ++ .../TPCInterpolationSpec.h | 49 ++++ .../TPCResidualWriterSpec.h | 53 ++++ .../TrackInterpolationReaderSpec.h | 26 ++ .../TrackInterpolationWorkflow.h | 28 ++ .../src/TPCInterpolationSpec.cxx | 258 ++++++++++++++++++ .../src/TPCResidualWriterSpec.cxx | 97 +++++++ .../src/TrackInterpolationReaderSpec.cxx | 21 ++ .../src/TrackInterpolationWorkflow.cxx | 61 +++++ .../src/tpc-interpolation-workflow.cxx | 42 +++ 11 files changed, 664 insertions(+), 1 deletion(-) create mode 100644 Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/CMakeLists.txt create mode 100644 Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCInterpolationSpec.h create mode 100644 Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualWriterSpec.h create mode 100644 Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TrackInterpolationReaderSpec.h create mode 100644 Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TrackInterpolationWorkflow.h create mode 100644 Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx create mode 100644 Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualWriterSpec.cxx create mode 100644 Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TrackInterpolationReaderSpec.cxx create mode 100644 Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TrackInterpolationWorkflow.cxx create mode 100644 Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-interpolation-workflow.cxx diff --git a/Detectors/GlobalTrackingWorkflow/CMakeLists.txt b/Detectors/GlobalTrackingWorkflow/CMakeLists.txt index 609896a533958..729d108fec904 100644 --- a/Detectors/GlobalTrackingWorkflow/CMakeLists.txt +++ b/Detectors/GlobalTrackingWorkflow/CMakeLists.txt @@ -22,4 +22,4 @@ o2_add_executable(match-workflow PUBLIC_LINK_LIBRARIES O2::GlobalTrackingWorkflow ) add_subdirectory(tofworkflow) - +add_subdirectory(tpcinterpolationworkflow) diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/CMakeLists.txt b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/CMakeLists.txt new file mode 100644 index 0000000000000..d602b72f92c7b --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/CMakeLists.txt @@ -0,0 +1,28 @@ +# Copyright CERN and copyright holders of ALICE O2. This software is distributed +# under the terms of the GNU General Public License v3 (GPL Version 3), copied +# verbatim in the file "COPYING". +# +# See http://alice-o2.web.cern.ch/license for full licensing information. +# +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization or +# submit itself to any jurisdiction. + +#TODO Does the O2::GlobalTracking library need to be linked? +o2_add_library(TPCInterpolationWorkflow + SOURCES src/TPCInterpolationSpec.cxx + src/TPCResidualWriterSpec.cxx + src/TrackInterpolationReaderSpec.cxx + src/TrackInterpolationWorkflow.cxx + PUBLIC_LINK_LIBRARIES O2::GlobalTracking + O2::ITSWorkflow + O2::SpacePoints + O2::GlobalTrackingWorkflow + O2::TOFWorkflow + O2::Framework + ) + +o2_add_executable(scdcalib-interpolation-workflow + COMPONENT_NAME tpc + SOURCES src/tpc-interpolation-workflow.cxx + PUBLIC_LINK_LIBRARIES O2::TPCInterpolationWorkflow) diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCInterpolationSpec.h b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCInterpolationSpec.h new file mode 100644 index 0000000000000..070686d8ed2b3 --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCInterpolationSpec.h @@ -0,0 +1,49 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_TPC_INTERPOLATION_SPEC_H +#define O2_TPC_INTERPOLATION_SPEC_H + +/// @file TPCInterpolationSpec.h + +#include "DataFormatsTPC/Constants.h" +#include "SpacePoints/TrackInterpolation.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" + +using namespace o2::framework; + +namespace o2 +{ +namespace tpc +{ +class TPCInterpolationDPL : public Task +{ + public: + TPCInterpolationDPL(bool useMC, const std::vector& tpcClusLanes) : mUseMC(useMC), mTPCClusLanes(tpcClusLanes) {} + ~TPCInterpolationDPL() override = default; + void init(InitContext& ic) final; + void run(ProcessingContext& pc) final; + + private: + o2::tpc::TrackInterpolation mInterpolation; // track interpolation engine + std::vector mTPCClusLanes; + std::array, o2::tpc::Constants::MAXSECTOR> mBufferedTPCClusters; + + bool mUseMC{false}; ///< MC flag +}; + +/// create a processor spec +framework::DataProcessorSpec getTPCInterpolationSpec(bool useMC, const std::vector& tpcClusLanes); + +} // namespace tpc +} // namespace o2 + +#endif diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualWriterSpec.h b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualWriterSpec.h new file mode 100644 index 0000000000000..1f9d7beb81b4f --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCResidualWriterSpec.h @@ -0,0 +1,53 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_TPC_RESIDUAL_WRITER_H +#define O2_TPC_RESIDUAL_WRITER_H + +/// @file TPCResidualWriterSpec.h + +#include "TFile.h" +#include "TTree.h" + +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" +#include + +namespace o2 +{ +namespace tpc +{ + +class ResidualWriterTPC : public o2::framework::Task +{ + public: + ResidualWriterTPC(bool useMC = false) : mUseMC(useMC) {} + ~ResidualWriterTPC() override = default; + void init(o2::framework::InitContext& ic) final; + void run(o2::framework::ProcessingContext& pc) final; + void endOfStream(o2::framework::EndOfStreamContext& ec) final; + + private: + bool mUseMC = false; ///< MC flag + std::string mOutFileName = "o2residuals_tpc.root"; ///< name of output file + std::string mTreeName = "residualsTPC"; ///< name of tree containing output + std::string mOutTracksBranchName = "tracks"; ///< name of branch containing output used tracks + std::string mOutResidualsBranchName = "residuals"; ///< name of branch containing output used residuals + std::unique_ptr mFile = nullptr; + std::unique_ptr mTree = nullptr; +}; + +/// create a processor spec +framework::DataProcessorSpec getTPCResidualWriterSpec(bool useMC); + +} // namespace tpc +} // namespace o2 + +#endif diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TrackInterpolationReaderSpec.h b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TrackInterpolationReaderSpec.h new file mode 100644 index 0000000000000..01896f687b6bf --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TrackInterpolationReaderSpec.h @@ -0,0 +1,26 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_TPC_INTERPOLATION_READER_H +#define O2_TPC_INTERPOLATION_READER_H + +/// @file TrackInterpolationReaderSpec.h + +namespace o2 +{ +namespace tpc +{ + +// TODO add specification for reading TPC cluster residuals and reference tracks + +} // namespace tpc +} // namespace o2 + +#endif diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TrackInterpolationWorkflow.h b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TrackInterpolationWorkflow.h new file mode 100644 index 0000000000000..2a611e7189d84 --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TrackInterpolationWorkflow.h @@ -0,0 +1,28 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_TPC_INTERPOLATION_WORKFLOW_H +#define O2_TPC_INTERPOLATION_WORKFLOW_H + +/// @file TrackInterpolationWorkflow.h + +#include "Framework/WorkflowSpec.h" + +namespace o2 +{ +namespace tpc +{ + +framework::WorkflowSpec getTPCInterpolationWorkflow(bool useMC); + +} // namespace tpc +} // namespace o2 + +#endif diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx new file mode 100644 index 0000000000000..852e60d5b6ef5 --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx @@ -0,0 +1,258 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @file TPCInterpolationSpec.cxx + +#include + +#include "Framework/ControlService.h" +#include "DataFormatsITS/TrackITS.h" +#include "ReconstructionDataFormats/TrackTPCITS.h" +#include "DataFormatsTPC/TrackTPC.h" +#include "DataFormatsTPC/ClusterNative.h" +#include "DataFormatsTPC/ClusterNativeHelper.h" +#include "DataFormatsTPC/TPCSectorHeader.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" +#include "TPCInterpolationWorkflow/TPCInterpolationSpec.h" + +using namespace o2::framework; + +namespace o2 +{ +namespace tpc +{ + +void TPCInterpolationDPL::init(InitContext& ic) +{ + //-------- init geometry and field --------// + o2::base::GeometryManager::loadGeometry(); + o2::base::Propagator::initFieldFromGRP("o2sim_grp.root"); + + mInterpolation.init(); +} + +void TPCInterpolationDPL::run(ProcessingContext& pc) +{ + + const auto tracksITS = pc.inputs().get>("trackITS"); + const auto tracksTPC = pc.inputs().get>("trackTPC"); + const auto tracksITSTPC = pc.inputs().get>("match"); + const auto tracksTPCClRefs = pc.inputs().get>("trackTPCClRefs"); + const auto trackMatchesTOF = pc.inputs().get>("matchTOF"); // FIXME missing reader + const auto clustersTOFInp = pc.inputs().get>("clustersTOF"); // FIXME o2::tof::Cluster is not messageable type which is required to create span + // make copy of TOF clusters... Is it needed? + std::vector clustersTOFCopy = clustersTOFInp; + auto clustersTOF = gsl::make_span(clustersTOFInp); + + // TPC Cluster loading part is copied from TPCITSMatchingSpec.cxx + //---------------------------->> TPC Clusters loading >>------------------------------------------ + int operation = 0; + uint64_t activeSectors = 0; + std::bitset validSectors = 0; + std::map datarefs; + for (auto const& lane : mTPCClusLanes) { + std::string inputLabel = "clusTPC" + std::to_string(lane); + LOG(INFO) << "Reading lane " << lane << " " << inputLabel; + auto ref = pc.inputs().get(inputLabel); + auto const* sectorHeader = DataRefUtils::getHeader(ref); + if (sectorHeader == nullptr) { + // FIXME: think about error policy + LOG(ERROR) << "sector header missing on header stack"; + return; + } + const int& sector = sectorHeader->sector; + // check the current operation, this is used to either signal eod or noop + // FIXME: the noop is not needed any more once the lane configuration with one + // channel per sector is used + if (sector < 0) { + if (operation < 0 && operation != sector) { + // we expect the same operation on all inputs + LOG(ERROR) << "inconsistent lane operation, got " << sector << ", expecting " << operation; + } else if (operation == 0) { + // store the operation + operation = sector; + } + continue; + } + if (validSectors.test(sector)) { + // have already data for this sector, this should not happen in the current + // sequential implementation, for parallel path merged at the tracker stage + // multiple buffers need to be handled + throw std::runtime_error("can only have one data set per sector"); + } + activeSectors |= sectorHeader->activeSectors; + validSectors.set(sector); + datarefs[sector] = ref; + } + + if (operation == -1) { + // EOD is transmitted in the sectorHeader with sector number equal to -1 + LOG(WARNING) << "operation = " << operation; + } + + auto printInputLog = [&validSectors, &activeSectors](auto& r, const char* comment, auto& s) { + LOG(INFO) << comment << " " << *(r.spec) << ", size " << DataRefUtils::getPayloadSize(r) // + << " for sector " << s // + << std::endl // + << " input status: " << validSectors // + << std::endl // + << " active sectors: " << std::bitset(activeSectors); + }; + + if (activeSectors == 0 || (activeSectors & validSectors.to_ulong()) != activeSectors) { + // not all sectors available + // Since we expect complete input, this should not happen (why does the bufferization considered for TPC CA tracker? Ask Matthias) + throw std::runtime_error("Did not receive TPC clusters data for all sectors"); + /* + for (auto const& refentry : datarefs) { + auto& sector = refentry.first; + auto& ref = refentry.second; + auto payploadSize = DataRefUtils::getPayloadSize(ref); + mBufferedTPCClusters[sector].resize(payploadSize); + std::copy(ref.payload, ref.payload + payploadSize, mBufferedTPCClusters[sector].begin()); + + printInputLog(ref, "buffering", sector); + } + // not needed to send something, DPL will simply drop this timeslice, whenever the + // data for all sectors is available, the output is sent in that time slice + return; + */ + } + //------------------------------------------------------------------------------ + std::array, o2::tpc::Constants::MAXSECTOR> mcInputs; // DUMMY + std::array, o2::tpc::Constants::MAXSECTOR> clustersTPC; + auto sectorStatus = validSectors; + + for (auto const& refentry : datarefs) { + auto& sector = refentry.first; + auto& ref = refentry.second; + clustersTPC[sector] = gsl::span(ref.payload, DataRefUtils::getPayloadSize(ref)); + sectorStatus.reset(sector); + printInputLog(ref, "received", sector); + } + if (sectorStatus.any()) { + LOG(ERROR) << "Reading bufferized TPC clusters, this should not happen"; + // some of the inputs have been buffered + for (size_t sector = 0; sector < sectorStatus.size(); ++sector) { + if (sectorStatus.test(sector)) { + clustersTPC[sector] = gsl::span(&mBufferedTPCClusters[sector].front(), mBufferedTPCClusters[sector].size()); + } + } + } + + // Just print TPC clusters status + { + // make human readable information from the bitfield + std::string bitInfo; + auto nActiveBits = validSectors.count(); + if (((uint64_t)0x1 << nActiveBits) == validSectors.to_ulong() + 1) { + // sectors 0 to some upper bound are active + bitInfo = "0-" + std::to_string(nActiveBits - 1); + } else { + int rangeStart = -1; + int rangeEnd = -1; + for (size_t sector = 0; sector < validSectors.size(); sector++) { + if (validSectors.test(sector)) { + if (rangeStart < 0) { + if (rangeEnd >= 0) { + bitInfo += ","; + } + bitInfo += std::to_string(sector); + if (nActiveBits == 1) { + break; + } + rangeStart = sector; + } + rangeEnd = sector; + } else { + if (rangeStart >= 0 && rangeEnd > rangeStart) { + bitInfo += "-" + std::to_string(rangeEnd); + } + rangeStart = -1; + } + } + if (rangeStart >= 0 && rangeEnd > rangeStart) { + bitInfo += "-" + std::to_string(rangeEnd); + } + } + LOG(INFO) << "running matching for sector(s) " << bitInfo; + } + + o2::tpc::ClusterNativeAccess clusterIndex; + std::unique_ptr clusterBuffer; + o2::tpc::MCLabelContainer clusterMCBuffer; + memset(&clusterIndex, 0, sizeof(clusterIndex)); + o2::tpc::ClusterNativeHelper::Reader::fillIndex(clusterIndex, clusterBuffer, clusterMCBuffer, clustersTPC, mcInputs, [&validSectors](auto& index) { return validSectors.test(index); }); + //----------------------------<< TPC Clusters loading <<------------------------------------------ + + // pass input data to TrackInterpolation object + mInterpolation.setITSTracksInp(tracksITS); + mInterpolation.setTPCTracksInp(tracksTPC); + mInterpolation.setTPCTrackClusIdxInp(tracksTPCClRefs); + mInterpolation.setTPCClustersInp(&clusterIndex); + mInterpolation.setTOFMatchesInp(trackMatchesTOF); + mInterpolation.setITSTPCTrackMatchesInp(tracksITSTPC); + mInterpolation.setTOFClustersInp(clustersTOF); + + if (mUseMC) { + // possibly MC labels will be used to check filtering procedure performance before interpolation + // not yet implemented + } + + printf("TPC Interpolation Workflow initialized. Start processing...\n"); + + mInterpolation.process(); + + pc.outputs().snapshot(Output{"GLO", "TPCINT_TRK", 0, Lifetime::Timeframe}, mInterpolation.getReferenceTracks()); + pc.outputs().snapshot(Output{"GLO", "TPCINT_RES", 0, Lifetime::Timeframe}, mInterpolation.getClusterResiduals()); +} + +DataProcessorSpec getTPCInterpolationSpec(bool useMC, const std::vector& tpcClusLanes) +{ + std::vector inputs; + std::vector outputs; + + inputs.emplace_back("trackITS", "ITS", "TRACKS", 0, Lifetime::Timeframe); + inputs.emplace_back("trackTPC", "TPC", "TRACKS", 0, Lifetime::Timeframe); + inputs.emplace_back("trackTPCClRefs", "TPC", "CLUSREFS", 0, Lifetime::Timeframe); + + for (auto lane : tpcClusLanes) { + std::string clusBind = "clusTPC" + std::to_string(lane); + inputs.emplace_back(clusBind.c_str(), "TPC", "CLUSTERNATIVE", lane, Lifetime::Timeframe); + } + + inputs.emplace_back("match", "GLO", "TPCITS", 0, Lifetime::Timeframe); + inputs.emplace_back("matchTOF", "TOF", "MATCHINFOS", 0, Lifetime::Timeframe); + inputs.emplace_back("clustersTOF", "TOF", "CLUSTERS", 0, Lifetime::Timeframe); + + if (useMC) { + LOG(FATAL) << "MC usage must be disabled for this workflow, since it is not yet implemented"; + // are the MC inputs from ITS-TPC matching and TOF matching duplicates? if trackITSMCTR == matchTOFMCITS one of them should be removed + inputs.emplace_back("trackITSMCTR", "GLO", "TPCITS_ITSMC", 0, Lifetime::Timeframe); + inputs.emplace_back("trackTPCMCTR", "GLO", "TPCITS_TPCMC", 0, Lifetime::Timeframe); + inputs.emplace_back("matchTOFMC", o2::header::gDataOriginTOF, "MATCHTOFINFOSMC", 0, Lifetime::Timeframe); + inputs.emplace_back("matchTOFMCTPC", o2::header::gDataOriginTOF, "MATCHTPCINFOSMC", 0, Lifetime::Timeframe); + inputs.emplace_back("matchTOFMCITS", o2::header::gDataOriginTOF, "MATCHITSINFOSMC", 0, Lifetime::Timeframe); + } + + outputs.emplace_back("GLO", "TPCINT_TRK", 0, Lifetime::Timeframe); + outputs.emplace_back("GLO", "TPCINT_RES", 0, Lifetime::Timeframe); + + return DataProcessorSpec{ + "tpc-track-interpolation", + inputs, + outputs, + AlgorithmSpec{adaptFromTask(useMC, tpcClusLanes)}, + Options{}}; +} + +} // namespace tpc +} // namespace o2 diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualWriterSpec.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualWriterSpec.cxx new file mode 100644 index 0000000000000..0f80013767ada --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualWriterSpec.cxx @@ -0,0 +1,97 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @file TPCResidualWriterSpec.cxx + +#include + +#include "Framework/ControlService.h" +#include "Framework/ConfigParamRegistry.h" +#include "SpacePoints/TrackInterpolation.h" +#include "TPCInterpolationWorkflow/TPCResidualWriterSpec.h" +#include "CommonUtils/StringUtils.h" + +using namespace o2::framework; + +namespace o2 +{ +namespace tpc +{ + +template +TBranch* getOrMakeBranch(TTree* tree, const char* brname, T* ptr) +{ + if (auto br = tree->GetBranch(brname)) { + // return branch if it already exists + br->SetAddress(static_cast(ptr)); + return br; + } + // otherwise make it + return tree->Branch(brname, ptr); +} + +void ResidualWriterTPC::init(InitContext& ic) +{ + mOutFileName = ic.options().get("residuals-outfile"); + mFile = std::make_unique(mOutFileName.c_str(), "RECREATE"); + if (!mFile->IsOpen()) { + throw std::runtime_error(o2::utils::concat_string("Failed to open TPC SCD residuals output file ", mOutFileName)); + } + mTree = std::make_unique(mTreeName.c_str(), "Tree of TPC residuals and reference tracks"); +} + +void ResidualWriterTPC::run(ProcessingContext& pc) +{ + + auto tracks = std::move(pc.inputs().get>("tracks")); + auto residuals = std::move(pc.inputs().get>("residuals")); + + auto tracksPtr = &tracks; + auto residualsPtr = &residuals; + + LOG(INFO) << "ResidualWriterTPC pulled " << tracks.size() << " reference tracks and " << residuals.size() << " TPC cluster residuals"; + + getOrMakeBranch(mTree.get(), mOutTracksBranchName.c_str(), &tracksPtr); + getOrMakeBranch(mTree.get(), mOutResidualsBranchName.c_str(), &residualsPtr); + + if (mUseMC) { + // TODO + } + + mTree->Fill(); +} + +void ResidualWriterTPC::endOfStream(EndOfStreamContext& ec) +{ + LOG(INFO) << "Finalizing TPC SCD interpolation residuals writing"; + mTree->Write(); + mTree.release()->Delete(); + mFile->Close(); +} + +DataProcessorSpec getTPCResidualWriterSpec(bool useMC) +{ + std::vector inputs; + inputs.emplace_back("tracks", "GLO", "TPCINT_TRK", 0, Lifetime::Timeframe); + inputs.emplace_back("residuals", "GLO", "TPCINT_RES", 0, Lifetime::Timeframe); + if (useMC) { + // TODO + } + return DataProcessorSpec{ + "tpc-residuals-writer", + inputs, + Outputs{}, + AlgorithmSpec{adaptFromTask(useMC)}, + Options{ + {"residuals-outfile", VariantType::String, "o2residuals_tpc.root", {"Name of the output file"}}}}; +} + +} // namespace tpc +} // namespace o2 diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TrackInterpolationReaderSpec.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TrackInterpolationReaderSpec.cxx new file mode 100644 index 0000000000000..c2e41024be2cf --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TrackInterpolationReaderSpec.cxx @@ -0,0 +1,21 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @file TrackInterpolationReaderSpec.cxx + +#include "TPCInterpolationWorkflow/TrackInterpolationReaderSpec.h" + +namespace o2 +{ +namespace tpc +{ + +} // namespace tpc +} // namespace o2 diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TrackInterpolationWorkflow.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TrackInterpolationWorkflow.cxx new file mode 100644 index 0000000000000..a326e26ee7352 --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TrackInterpolationWorkflow.cxx @@ -0,0 +1,61 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @file TrackInterpolationWorkflow.cxx + +#include + +#include "ITSWorkflow/TrackReaderSpec.h" +#include "TPCWorkflow/TrackReaderSpec.h" +#include "TPCWorkflow/PublisherSpec.h" +#include "GlobalTrackingWorkflow/TrackTPCITSReaderSpec.h" +#include "TOFWorkflow/ClusterReaderSpec.h" +#include "TOFWorkflow/TOFMatchedReaderSpec.h" +#include "Algorithm/RangeTokenizer.h" +#include "TPCInterpolationWorkflow/TPCResidualWriterSpec.h" +#include "TPCInterpolationWorkflow/TPCInterpolationSpec.h" +#include "TPCInterpolationWorkflow/TrackInterpolationWorkflow.h" + +namespace o2 +{ +namespace tpc +{ + +framework::WorkflowSpec getTPCInterpolationWorkflow(bool useMC) +{ + framework::WorkflowSpec specs; + + specs.emplace_back(o2::its::getITSTrackReaderSpec(useMC)); + specs.emplace_back(o2::tpc::getTPCTrackReaderSpec(useMC)); + + std::vector tpcClusSectors = o2::RangeTokenizer::tokenize("0-35"); + std::vector tpcClusLanes = tpcClusSectors; + specs.emplace_back(o2::tpc::getPublisherSpec(o2::tpc::PublisherConf{ + "tpc-native-cluster-reader", + "tpcrec", + {"clusterbranch", "TPCClusterNative", "Branch with TPC native clusters"}, + {"clustermcbranch", "TPCClusterNativeMCTruth", "MC label branch"}, + OutputSpec{"TPC", "CLUSTERNATIVE"}, + OutputSpec{"TPC", "CLNATIVEMCLBL"}, + tpcClusSectors, + tpcClusLanes}, + useMC)); + + specs.emplace_back(o2::globaltracking::getTrackTPCITSReaderSpec(useMC)); + specs.emplace_back(o2::tof::getClusterReaderSpec(useMC)); + specs.emplace_back(o2::tof::getTOFMatchedReaderSpec(useMC)); + specs.emplace_back(o2::tpc::getTPCInterpolationSpec(useMC, tpcClusLanes)); + specs.emplace_back(o2::tpc::getTPCResidualWriterSpec(useMC)); + + return specs; +} + +} // namespace tpc +} // namespace o2 diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-interpolation-workflow.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-interpolation-workflow.cxx new file mode 100644 index 0000000000000..eb1b91357f146 --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-interpolation-workflow.cxx @@ -0,0 +1,42 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "TPCInterpolationWorkflow/TrackInterpolationWorkflow.h" +#include "CommonUtils/ConfigurableParam.h" + +using namespace o2::framework; + +// ------------------------------------------------------------------ + +// we need to add workflow options before including Framework/runDataProcessing +void customize(std::vector& workflowOptions) +{ + // option allowing to set parameters + workflowOptions.push_back(ConfigParamSpec{ + "disable-mc", o2::framework::VariantType::Bool, false, {"disable MC propagation even if available"}}); + + std::string keyvaluehelp("Semicolon separated key=value strings ..."); + workflowOptions.push_back(ConfigParamSpec{"configKeyValues", VariantType::String, "", {keyvaluehelp}}); +} + +// ------------------------------------------------------------------ + +#include "Framework/runDataProcessing.h" + +WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) +{ + // Update the (declared) parameters if changed from the command line + o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); + // write the configuration used for the workflow + o2::conf::ConfigurableParam::writeINI("o2tpcinterpolation-workflow_configuration.ini"); + + auto useMC = !configcontext.options().get("disable-mc"); + return std::move(o2::tpc::getTPCInterpolationWorkflow(useMC)); +}