From bb401f7106171a2d0e1cdf128a75a7b3a9348085 Mon Sep 17 00:00:00 2001 From: Matteo Concas Date: Wed, 21 Jul 2021 23:28:47 +0200 Subject: [PATCH 01/10] Add smoother files --- Detectors/ITSMFT/ITS/tracking/CMakeLists.txt | 1 + .../tracking/include/ITStracking/Smoother.h | 68 ++++++ .../ITSMFT/ITS/tracking/src/Smoother.cxx | 231 ++++++++++++++++++ macro/run_trac_ca_its.C | 15 +- 4 files changed, 311 insertions(+), 4 deletions(-) create mode 100644 Detectors/ITSMFT/ITS/tracking/include/ITStracking/Smoother.h create mode 100644 Detectors/ITSMFT/ITS/tracking/src/Smoother.cxx diff --git a/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt b/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt index 8e258fdf8dab6..bf7dcd8bf1364 100644 --- a/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt @@ -26,6 +26,7 @@ o2_add_library(ITStracking src/ClusterLines.cxx src/Vertexer.cxx src/VertexerTraits.cxx + src/Smoother.cxx PUBLIC_LINK_LIBRARIES O2::GPUCommon Microsoft.GSL::GSL O2::CommonConstants diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Smoother.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Smoother.h new file mode 100644 index 0000000000000..6d2a512bae19e --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Smoother.h @@ -0,0 +1,68 @@ +// 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 Smoother.h +/// \brief Class to handle Kalman smoothing for ITS tracking. +/// Its instance stores the state of the track to the level we want to smooth to avoid multiple re-propagations when testing different clusters. +/// + +#include "ReconstructionDataFormats/Track.h" +#include "DataFormatsITS/TrackITS.h" +#include "DetectorsBase/Propagator.h" +#include "ITStracking/ROframe.h" + +#if defined(CA_DEBUG) || defined(CA_STANDALONE_DEBUGGER) +#include "ITStracking/StandaloneDebugger.h" +#endif + +namespace o2 +{ +namespace its +{ + +template +class Smoother +{ + public: + Smoother(TrackITSExt& track, int layer, const ROframe& event, float bZ, o2::base::PropagatorF::MatCorrType corr); + ~Smoother(); + + bool isValidInit() const + { + return mInitStatus; + } + bool testCluster(const int clusterId, const ROframe& event); + bool getSmoothedTrack(); + float getChi2() const { return mBestChi2; } + float getLastChi2() const { return mLastChi2; } + + private: + float computeSmoothedPredictedChi2(const o2::track::TrackParCov& outwTrack, + const o2::track::TrackParCov& inwTrack, + const std::array& cls, + const std::array& clCov); + bool smoothTrack(); + + private: + int mLayerToSmooth; + float mBz; // Magnetic field along Z + bool mInitStatus; // State after the initialization + o2::base::PropagatorF::MatCorrType mCorr; // Type of correction to use + TrackITSExt mInwardsTrack; // outwards track: from innermost cluster to outermost + TrackITSExt mOutwardsTrack; // inwards track: from outermost cluster to innermost + float mBestChi2; // Best value of local smoothed chi2 + float mLastChi2 = 1e8; // Latest computed chi2 + +#if defined(CA_DEBUG) || defined(CA_STANDALONE_DEBUGGER) + StandaloneDebugger* mDebugger; +#endif +}; +} // namespace its +} // namespace o2 \ No newline at end of file diff --git a/Detectors/ITSMFT/ITS/tracking/src/Smoother.cxx b/Detectors/ITSMFT/ITS/tracking/src/Smoother.cxx new file mode 100644 index 0000000000000..6743ef2e865c1 --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/src/Smoother.cxx @@ -0,0 +1,231 @@ +// 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 "ITStracking/Smoother.h" + +namespace o2 +{ +namespace its +{ + +constexpr std::array getInverseSymm2D(const std::array& mat) +{ + const double det = mat[0] * mat[2] - mat[1] * mat[1]; + return std::array{mat[2] / det, -mat[1] / det, mat[0] / det}; +} + +// Smoother +template +Smoother::Smoother(TrackITSExt& track, int smoothingLayer, const ROframe& event, float bZ, o2::base::PropagatorF::MatCorrType corr) : mLayerToSmooth{smoothingLayer}, + mBz(bZ), + mCorr(corr) +{ +#if defined(CA_DEBUG) || defined(CA_STANDALONE_DEBUGGER) + mDebugger = new StandaloneDebugger("dbg_ITSSmootherCPU.root"); +#endif + + auto propInstance = o2::base::Propagator::Instance(); + const TrackingFrameInfo& originalTf = event.getTrackingFrameInfoOnLayer(mLayerToSmooth).at(track.getClusterIndex(mLayerToSmooth)); + + mOutwardsTrack = track; // This track will be propagated outwards inside the smoother! (as last step of fitting did inward propagation) + mInwardsTrack = {track.getParamOut(), // This track will be propagated inwards inside the smoother! + static_cast(mOutwardsTrack.getNumberOfClusters()), -999, static_cast(event.getROFrameId()), + mOutwardsTrack.getParamOut(), mOutwardsTrack.getClusterIndexes()}; + + mOutwardsTrack.resetCovariance(); + mOutwardsTrack.setChi2(0); + mInwardsTrack.resetCovariance(); + mInwardsTrack.setChi2(0); + + bool statusOutw{false}; + bool statusInw{false}; + + ////////////////////// + // Outward propagation + for (size_t iLayer{0}; iLayer < mLayerToSmooth; ++iLayer) { + if (mOutwardsTrack.getClusterIndex(iLayer) == constants::its::UnusedIndex) { // Shorter tracks + continue; + } + const TrackingFrameInfo& tF = event.getTrackingFrameInfoOnLayer(iLayer).at(mOutwardsTrack.getClusterIndex(iLayer)); + statusOutw = mOutwardsTrack.rotate(tF.alphaTrackingFrame); + statusOutw &= propInstance->propagateToX(mOutwardsTrack, + tF.xTrackingFrame, + mBz, + o2::base::PropagatorImpl::MAX_SIN_PHI, + o2::base::PropagatorImpl::MAX_STEP, + mCorr); + mOutwardsTrack.setChi2(mOutwardsTrack.getChi2() + mOutwardsTrack.getPredictedChi2(tF.positionTrackingFrame, tF.covarianceTrackingFrame)); + statusOutw &= mOutwardsTrack.o2::track::TrackParCov::update(tF.positionTrackingFrame, tF.covarianceTrackingFrame); + // LOG(INFO) << "Outwards loop on inwards track, layer: " << iLayer << " x: " << mOutwardsTrack.getX(); + } + + // Prediction on the previously outwards-propagated track is done on a copy, as the process seems to be not reversible + auto outwardsClone = mOutwardsTrack; + statusOutw = outwardsClone.rotate(originalTf.alphaTrackingFrame); + statusOutw &= propInstance->propagateToX(outwardsClone, + originalTf.xTrackingFrame, + mBz, + o2::base::PropagatorImpl::MAX_SIN_PHI, + o2::base::PropagatorImpl::MAX_STEP, + mCorr); + ///////////////////// + // Inward propagation + for (size_t iLayer{D - 1}; iLayer > mLayerToSmooth; --iLayer) { + if (mInwardsTrack.getClusterIndex(iLayer) == constants::its::UnusedIndex) { // Shorter tracks + continue; + } + const TrackingFrameInfo& tF = event.getTrackingFrameInfoOnLayer(iLayer).at(mInwardsTrack.getClusterIndex(iLayer)); + statusInw = mInwardsTrack.rotate(tF.alphaTrackingFrame); + statusInw &= propInstance->propagateToX(mInwardsTrack, + tF.xTrackingFrame, + mBz, + o2::base::PropagatorImpl::MAX_SIN_PHI, + o2::base::PropagatorImpl::MAX_STEP, + mCorr); + mInwardsTrack.setChi2(mInwardsTrack.getChi2() + mInwardsTrack.getPredictedChi2(tF.positionTrackingFrame, tF.covarianceTrackingFrame)); + statusInw &= mInwardsTrack.o2::track::TrackParCov::update(tF.positionTrackingFrame, tF.covarianceTrackingFrame); + // LOG(INFO) << "Inwards loop on outwards track, layer: " << iLayer << " x: " << mInwardsTrack.getX(); + } + + // Prediction on the previously inwards-propagated track is done on a copy, as the process seems to be not revesible + auto inwardsClone = mInwardsTrack; + statusInw = inwardsClone.rotate(originalTf.alphaTrackingFrame); + statusInw &= propInstance->propagateToX(inwardsClone, + originalTf.xTrackingFrame, + mBz, + o2::base::PropagatorImpl::MAX_SIN_PHI, + o2::base::PropagatorImpl::MAX_STEP, + mCorr); + // Compute weighted local chi2 + mInitStatus = statusInw && statusOutw; + if (mInitStatus) { + mBestChi2 = computeSmoothedPredictedChi2(inwardsClone, outwardsClone, originalTf.positionTrackingFrame, originalTf.covarianceTrackingFrame); + mLastChi2 = mBestChi2; + LOG(INFO) << "Smoothed chi2 on original cluster: " << mBestChi2; + } +} + +#if defined(CA_DEBUG) || defined(CA_STANDALONE_DEBUGGER) +template +Smoother::~Smoother() +{ + delete mDebugger; +} +#else +template +Smoother::~Smoother() = default; +#endif + +template +float Smoother::computeSmoothedPredictedChi2(const o2::track::TrackParCov& firstTrack, // outwards track: from innermost cluster to outermost + const o2::track::TrackParCov& secondTrack, // inwards track: from outermost cluster to innermost + const std::array& cls, + const std::array& clCov) +{ + // Tracks need to be already propagated, compute only chi2 + // Symmetric covariances assumed + + if (firstTrack.getX() != secondTrack.getX()) { + LOG(FATAL) << "Tracks need to be propagated to the same point! secondTrack.X=" << secondTrack.getX() << " firstTrack.X=" << firstTrack.getX(); + } + + std::array pp1 = {static_cast(firstTrack.getY()), static_cast(firstTrack.getZ())}; // P1: predicted Y,Z points + std::array pp2 = {static_cast(secondTrack.getY()), static_cast(secondTrack.getZ())}; // P2: predicted Y,Z points + + std::array c1 = {static_cast(firstTrack.getSigmaY2()), + static_cast(firstTrack.getSigmaZY()), + static_cast(firstTrack.getSigmaZ2())}; // Cov. track 1 + + std::array c2 = {static_cast(secondTrack.getSigmaY2()), + static_cast(secondTrack.getSigmaZY()), + static_cast(secondTrack.getSigmaZ2())}; // Cov. track 2 + + std::array w1 = getInverseSymm2D(c1); // weight matrices + std::array w2 = getInverseSymm2D(c2); + + std::array w1w2 = {w1[0] + w2[0], w1[1] + w2[1], w1[2] + w2[2]}; // (W1 + W2) + std::array C = getInverseSymm2D(w1w2); // C = (W1+W2)^-1 + + std::array w1pp1 = {w1[0] * pp1[0] + w1[1] * pp1[1], w1[1] * pp1[0] + w1[2] * pp1[1]}; // W1 * P1 + std::array w2pp2 = {w2[0] * pp2[0] + w2[1] * pp2[1], w2[1] * pp2[0] + w2[2] * pp2[1]}; // W2 * P2 + + double Y = C[0] * (w1pp1[0] + w2pp2[0]) + C[1] * (w1pp1[1] + w2pp2[1]); // Pp: weighted normalized combination of the predictions: + double Z = C[1] * (w1pp1[0] + w2pp2[0]) + C[2] * (w1pp1[1] + w2pp2[1]); // Pp = [(W1 * P1) + (W2 * P2)] / (W1 + W2) + + std::array delta = {Y - cls[0], Z - cls[1]}; // Δ = Pp - X, X: space point of cluster (Y,Z) + std::array CCp = {C[0] + static_cast(clCov[0]), C[1] + static_cast(clCov[1]), C[2] + static_cast(clCov[2])}; // Transformation of cluster covmat: CCp = C + Cov + std::array Wp = getInverseSymm2D(CCp); // Get weight matrix: Wp = CCp^-1 + + float chi2 = static_cast(delta[0] * (Wp[0] * delta[0] + Wp[1] * delta[1]) + delta[1] * (Wp[1] * delta[0] + Wp[2] * delta[1])); // chi2 = tΔ * (Wp * Δ) + + // #ifdef CA_DEBUG + LOG(INFO) << "Cluster_y: " << cls[0] << " Cluster_z: " << cls[1]; + LOG(INFO) << "\t\t- Covariance cluster: Y2: " << clCov[0] << " YZ: " << clCov[1] << " Z2: " << clCov[2]; + LOG(INFO) << "\t\t- Propagated t1_y: " << pp1[0] << " t1_z: " << pp1[1]; + LOG(INFO) << "\t\t- Propagated t2_y: " << pp2[0] << " t2_z: " << pp2[1]; + LOG(INFO) << "\t\t- Covariance t1: sY2: " << c1[0] << " sYZ: " << c1[1] << " sZ2: " << c1[2]; + LOG(INFO) << "\t\t- Covariance t2: sY2: " << c2[0] << " sYZ: " << c2[1] << " sZ2: " << c2[2]; + LOG(INFO) << "Smoother prediction Y: " << Y << " Z: " << Z; + LOG(INFO) << "\t\t- Delta_y: " << delta[0] << " Delta_z: " << delta[1]; + LOG(INFO) << "\t\t- Covariance Pr: Y2: " << C[0] << " YZ: " << C[1] << " Z2: " << C[2]; + LOG(INFO) << "\t\t- predicted chi2 t1: " << firstTrack.getPredictedChi2(cls, clCov); + LOG(INFO) << "\t\t- predicted chi2 t2: " << secondTrack.getPredictedChi2(cls, clCov); + // #endif + return chi2; +} + +template +bool Smoother::testCluster(const int clusterId, const ROframe& event) +{ + if (!mInitStatus) { + return false; + } + auto propInstance = o2::base::Propagator::Instance(); + const TrackingFrameInfo& testTf = event.getTrackingFrameInfoOnLayer(mLayerToSmooth).at(clusterId); + + bool statusOutw{false}; + bool statusInw{false}; + + // Prediction on the previously outwards-propagated track is done on a copy, as the process seems to be not reversible + auto outwardsClone = mOutwardsTrack; + statusOutw = outwardsClone.rotate(testTf.alphaTrackingFrame); + statusOutw &= propInstance->propagateToX(outwardsClone, + testTf.xTrackingFrame, + mBz, + o2::base::PropagatorImpl::MAX_SIN_PHI, + o2::base::PropagatorImpl::MAX_STEP, + mCorr); + + // Prediction on the previously inwards-propagated track is done on a copy, as the process seems to be not reversible + auto inwardsClone = mInwardsTrack; + statusInw = inwardsClone.rotate(testTf.alphaTrackingFrame); + statusInw &= propInstance->propagateToX(inwardsClone, + testTf.xTrackingFrame, + mBz, + o2::base::PropagatorImpl::MAX_SIN_PHI, + o2::base::PropagatorImpl::MAX_STEP, + mCorr); + bool testStatus = statusOutw && statusInw; + if (!(statusOutw && statusInw)) { + LOG(WARNING) << "Failed propagation in smoother!"; + return false; + } + + // Compute weighted local chi2 + mLastChi2 = computeSmoothedPredictedChi2(inwardsClone, outwardsClone, testTf.positionTrackingFrame, testTf.covarianceTrackingFrame); + LOG(INFO) << "Smoothed chi2 on tested cluster: " << mLastChi2; + + return true; +} + +template class Smoother<7>; + +} // namespace its +} // namespace o2 \ No newline at end of file diff --git a/macro/run_trac_ca_its.C b/macro/run_trac_ca_its.C index 129adba7bfaf5..621cdacc673fa 100644 --- a/macro/run_trac_ca_its.C +++ b/macro/run_trac_ca_its.C @@ -144,7 +144,8 @@ void run_trac_ca_its(bool cosmics = false, o2::itsmft::TopologyDictionary dict; if (dictfile.empty()) { - dictfile = o2::base::NameConf::getAlpideClusterDictionaryFileName(o2::detectors::DetID::ITS, "", "bin"); + dictfile = o2::base::NameConf::getDictionaryFileName(o2::detectors::DetID::ITS, "", "bin"); + std::cout << dictfile; } std::ifstream file(dictfile.c_str()); if (file.good()) { @@ -205,7 +206,7 @@ void run_trac_ca_its(bool cosmics = false, memParams[0].CellsMemoryCoefficients[iLayer] = 0.001f; } } else { - // PbPb tracking params + // Comment for pp // ---- trackParams.resize(3); memParams.resize(3); @@ -248,7 +249,11 @@ void run_trac_ca_its(bool cosmics = false, pattIt = patt.begin(); int rofId{0}; for (auto& rof : *rofs) { - + counterEvent++; + // if (counterEvent != 111) { + // continue; + // } + std::cout << "ROF: " << counterEvent << std::endl; auto start = std::chrono::high_resolution_clock::now(); auto it = pattIt; o2::its::ioutils::loadROFrameData(rof, event, clSpan, pattIt, dict, labels); @@ -270,7 +275,9 @@ void run_trac_ca_its(bool cosmics = false, if (!vertITS.empty()) { // Using only the first vertex in the list std::cout << " - Reconstructed vertex: x = " << vertITS[0].getX() << " y = " << vertITS[0].getY() << " x = " << vertITS[0].getZ() << std::endl; - event.addPrimaryVertex(vertITS[0].getX(), vertITS[0].getY(), vertITS[0].getZ()); + // LOG(FATAL) << event.getPrimaryVertex(0).x << " " << event.getPrimaryVertex(0).y << " " << event.getPrimaryVertex(0).z; + // LOG(FATAL) << " --- " << event.getPrimaryVerticesNum(); + // event.addPrimaryVertex(vertITS[0].getX(), vertITS[0].getY(), vertITS[0].getZ()); } else { std::cout << " - Vertex not reconstructed, tracking skipped" << std::endl; } From b8292cb7ccc2d6319d3e2bcd513f9cc6fd34d942 Mon Sep 17 00:00:00 2001 From: Matteo Concas Date: Thu, 29 Jul 2021 16:59:51 +0200 Subject: [PATCH 02/10] Add smoother header in tracker --- .../ITSMFT/ITS/tracking/include/ITStracking/Smoother.h | 2 +- Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h | 6 ++++++ Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx | 4 ++++ 3 files changed, 11 insertions(+), 1 deletion(-) diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Smoother.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Smoother.h index 6d2a512bae19e..6d075cead0421 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Smoother.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Smoother.h @@ -51,7 +51,7 @@ class Smoother bool smoothTrack(); private: - int mLayerToSmooth; + int mLayerToSmooth; // Layer to compute smoothing optimization float mBz; // Magnetic field along Z bool mInitStatus; // State after the initialization o2::base::PropagatorF::MatCorrType mCorr; // Type of correction to use diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h index 47fdadb44a0cb..3e3d80f5cb2a7 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h @@ -70,6 +70,11 @@ class Tracker float getBz() const; void clustersToTracks(std::function = [](std::string s) { std::cout << s << std::endl; }); + void setSmoothing(bool v) { mApplySmoothing = v; } + bool getSmoothing() const { return mApplySmoothing; } + + std::vector& getTracks(); + auto& getTrackLabels() { return mTrackLabels; } void setROFrame(std::uint32_t f) { mROFrame = f; } std::uint32_t getROFrame() const { return mROFrame; } @@ -104,6 +109,7 @@ class Tracker std::vector mTrkParams; bool mCUDA = false; + bool mApplySmoothing = false; o2::base::PropagatorImpl::MatCorrType mCorrType = o2::base::PropagatorImpl::MatCorrType::USEMatCorrLUT; float mBz = 5.f; std::uint32_t mROFrame = 0; diff --git a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx index 5fdbfb088dcbd..d5dbe3a2c353a 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx @@ -18,6 +18,7 @@ #include "ITStracking/Cell.h" #include "ITStracking/Constants.h" #include "ITStracking/IndexTableUtils.h" +#include "ITStracking/Smoother.h" #include "ITStracking/Tracklet.h" #include "ITStracking/TrackerTraits.h" #include "ITStracking/TrackerTraitsCPU.h" @@ -307,6 +308,9 @@ void Tracker::findTracks() tracks.emplace_back(temporaryTrack); } + if (mApplySmoothing) { + // Smoothing tracks + } std::sort(tracks.begin(), tracks.end(), [](TrackITSExt& track1, TrackITSExt& track2) { return track1.isBetter(track2, 1.e6f); }); From d9782642213175e20e4c6ff4165fd698057f4216 Mon Sep 17 00:00:00 2001 From: Matteo Concas Date: Tue, 10 Aug 2021 19:29:08 +0200 Subject: [PATCH 03/10] Add macro to test MC clusters --- Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt | 7 +++++++ .../ITSMFT/ITS/macros/test/CheckMCTracks.C | 17 +++++++++++++++++ 2 files changed, 24 insertions(+) create mode 100644 Detectors/ITSMFT/ITS/macros/test/CheckMCTracks.C diff --git a/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt b/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt index befc776289e25..0e9b44d468a0a 100644 --- a/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt @@ -70,6 +70,13 @@ o2_add_test_root_macro(RunGPUTracking.C PUBLIC_LINK_LIBRARIES O2::GPUTracking LABELS its) +o2_add_test_root_macro(CheckMCTracks.C + PUBLIC_LINK_LIBRARIES O2::SimulationDataFormat + O2::ITSBase + O2::DataFormatsITS + O2::DataFormatsITSMFT + LABELS its) + o2_add_test_root_macro(DisplayTrack.C PUBLIC_LINK_LIBRARIES O2::ITSBase O2::DataFormatsITSMFT diff --git a/Detectors/ITSMFT/ITS/macros/test/CheckMCTracks.C b/Detectors/ITSMFT/ITS/macros/test/CheckMCTracks.C new file mode 100644 index 0000000000000..51ad9dcaa9bd6 --- /dev/null +++ b/Detectors/ITSMFT/ITS/macros/test/CheckMCTracks.C @@ -0,0 +1,17 @@ +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#endif + +void CheckMCTrack() +{ + // Reconstructed tracks + TFile* file1 = TFile::Open(tracfile.data()); + TTree* recTree = (TTree*)gFile->Get("o2sim"); + std::vector* recArr = nullptr; + recTree->SetBranchAddress("ITSTrack", &recArr); + // Track MC labels + std::vector* trkLabArr = nullptr; + recTree->SetBranchAddress("ITSTrackMCTruth", &trkLabArr); +} \ No newline at end of file From cd731c0e031f1eee9e00ad55ed28bbce3b675ff5 Mon Sep 17 00:00:00 2001 From: Matteo Concas Date: Wed, 11 Aug 2021 17:48:47 +0200 Subject: [PATCH 04/10] Add macro to check fake clusters on reco tracks --- .../ITS/include/DataFormatsITS/TrackITS.h | 1 + .../Detectors/ITSMFT/ITS/src/TrackITS.cxx | 11 +++ .../ITSMFT/ITS/macros/test/CMakeLists.txt | 2 +- .../ITSMFT/ITS/macros/test/CheckMCTracks.C | 17 ----- .../macros/test/CheckTrackerFakeClusters.C | 69 +++++++++++++++++++ .../tracking/include/ITStracking/Tracker.h | 1 - 6 files changed, 82 insertions(+), 19 deletions(-) delete mode 100644 Detectors/ITSMFT/ITS/macros/test/CheckMCTracks.C create mode 100644 Detectors/ITSMFT/ITS/macros/test/CheckTrackerFakeClusters.C diff --git a/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h b/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h index a4b1510e0a289..6968aa8ade4c4 100644 --- a/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h +++ b/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h @@ -98,6 +98,7 @@ class TrackITS : public o2::track::TrackParCov int getPattern() const { return mPattern; } bool hasHitOnLayer(int i) { return mPattern & (0x1 << i); } bool isFakeOnLayer(int i) { return !(mPattern & (0x1 << (16 + i))); } + int getNFakeClusters(); void setNextROFbit(bool toggle = true) { setUserField((getUserField() & ~kNextROF) | (-toggle & kNextROF)); } bool hasHitInNextROF() const { return getUserField() & kNextROF; } diff --git a/DataFormats/Detectors/ITSMFT/ITS/src/TrackITS.cxx b/DataFormats/Detectors/ITSMFT/ITS/src/TrackITS.cxx index 9b2d3563b45a9..6f0fb6256daf2 100644 --- a/DataFormats/Detectors/ITSMFT/ITS/src/TrackITS.cxx +++ b/DataFormats/Detectors/ITSMFT/ITS/src/TrackITS.cxx @@ -106,3 +106,14 @@ Bool_t TrackITS::isBetter(const TrackITS& best, Float_t maxChi2) const } return kFALSE; } + +int TrackITS::getNFakeClusters() +{ + int nFake{0}; + for (int iCl{0}; iCl < getNClusters(); ++iCl) { + if (hasHitOnLayer(iCl) && isFakeOnLayer(iCl)) { + ++nFake; + } + } + return nFake; +} \ No newline at end of file diff --git a/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt b/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt index 0e9b44d468a0a..dcd6026bd9094 100644 --- a/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt @@ -70,7 +70,7 @@ o2_add_test_root_macro(RunGPUTracking.C PUBLIC_LINK_LIBRARIES O2::GPUTracking LABELS its) -o2_add_test_root_macro(CheckMCTracks.C +o2_add_test_root_macro(CheckTrackerFakeClusters.C PUBLIC_LINK_LIBRARIES O2::SimulationDataFormat O2::ITSBase O2::DataFormatsITS diff --git a/Detectors/ITSMFT/ITS/macros/test/CheckMCTracks.C b/Detectors/ITSMFT/ITS/macros/test/CheckMCTracks.C deleted file mode 100644 index 51ad9dcaa9bd6..0000000000000 --- a/Detectors/ITSMFT/ITS/macros/test/CheckMCTracks.C +++ /dev/null @@ -1,17 +0,0 @@ -#if !defined(__CLING__) || defined(__ROOTCLING__) -#include -#include -#include -#endif - -void CheckMCTrack() -{ - // Reconstructed tracks - TFile* file1 = TFile::Open(tracfile.data()); - TTree* recTree = (TTree*)gFile->Get("o2sim"); - std::vector* recArr = nullptr; - recTree->SetBranchAddress("ITSTrack", &recArr); - // Track MC labels - std::vector* trkLabArr = nullptr; - recTree->SetBranchAddress("ITSTrackMCTruth", &trkLabArr); -} \ No newline at end of file diff --git a/Detectors/ITSMFT/ITS/macros/test/CheckTrackerFakeClusters.C b/Detectors/ITSMFT/ITS/macros/test/CheckTrackerFakeClusters.C new file mode 100644 index 0000000000000..0a5b62c1212d6 --- /dev/null +++ b/Detectors/ITSMFT/ITS/macros/test/CheckTrackerFakeClusters.C @@ -0,0 +1,69 @@ +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "DataFormatsITS/TrackITS.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#endif + +void CheckMCTracks(const std::string tracfile = "o2trac_its.root") +{ + TFile* file1 = TFile::Open(tracfile.data()); + TTree* recTree = (TTree*)file1->Get("o2sim"); + std::vector* recArr = nullptr; + recTree->SetBranchAddress("ITSTrack", &recArr); + std::vector* trkLabArr = nullptr; + recTree->SetBranchAddress("ITSTrackMCTruth", &trkLabArr); + std::vector histLength; + std::vector histLength1Fake; + histLength.resize(4); + histLength1Fake.resize(4); + for (int iH{4}; iH < 8; ++iH) { + histLength[iH - 4] = new TH1I(Form("trk_len_%d", iH), Form("trk_len = %d", iH), 7, -.5, 6.5); + } + for (int iH{4}; iH < 8; ++iH) { + histLength1Fake[iH - 4] = new TH1I(Form("trk_len_%d_1f", iH), Form("trk_len = %d, 1 Fake", iH), 7, -.5, 6.5); + } + auto nt = new TNtuple("TrackITSInfo", "Info about ITS tracks", "len:l0:l1:l2:l3:l4:l5:6"); + std::array labels; + for (unsigned int iEntry{0}; iEntry < recTree->GetEntriesFast(); ++iEntry) { + recTree->GetEntry(iEntry); + for (unsigned int iTrack{0}; iTrack < recArr->size(); ++iTrack) { + auto& track = (*recArr)[iTrack]; + auto len = (*recArr)[iTrack].getNClusters(); + int count = 0; + for (int iLabel{0}; iLabel < 7; ++iLabel) { + if (track.hasHitOnLayer(iLabel)) { + count++; + if (track.isFakeOnLayer(iLabel)) { + histLength[len - 4]->Fill(iLabel); + if (track.getNFakeClusters() == 1) { + histLength1Fake[len - 4]->Fill(iLabel); + } + } + } + } + if (count != len) { + std::cout << "mismatch!\n"; + } + } + } + auto canvas = new TCanvas("c1", "c1", 2400, 1600); + canvas->Divide(4, 2); + for (int iH{0}; iH < 4; ++iH) { + canvas->cd(iH + 1); + histLength[iH]->Draw(); + } + for (int iH{0}; iH < 4; ++iH) { + canvas->cd(iH + 5); + histLength1Fake[iH]->Draw(); + } + canvas->SaveAs("fakeClusters.png", "recreate"); +} \ No newline at end of file diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h index 3e3d80f5cb2a7..1c4bb4f542c87 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h @@ -74,7 +74,6 @@ class Tracker bool getSmoothing() const { return mApplySmoothing; } std::vector& getTracks(); - auto& getTrackLabels() { return mTrackLabels; } void setROFrame(std::uint32_t f) { mROFrame = f; } std::uint32_t getROFrame() const { return mROFrame; } From d8d94bdd267b9f035c06c9dc33c279899fb34c0b Mon Sep 17 00:00:00 2001 From: Matteo Concas Date: Tue, 14 Sep 2021 11:07:33 +0200 Subject: [PATCH 05/10] Update run_trac_ca_its.C --- macro/run_trac_ca_its.C | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/macro/run_trac_ca_its.C b/macro/run_trac_ca_its.C index 621cdacc673fa..129adba7bfaf5 100644 --- a/macro/run_trac_ca_its.C +++ b/macro/run_trac_ca_its.C @@ -144,8 +144,7 @@ void run_trac_ca_its(bool cosmics = false, o2::itsmft::TopologyDictionary dict; if (dictfile.empty()) { - dictfile = o2::base::NameConf::getDictionaryFileName(o2::detectors::DetID::ITS, "", "bin"); - std::cout << dictfile; + dictfile = o2::base::NameConf::getAlpideClusterDictionaryFileName(o2::detectors::DetID::ITS, "", "bin"); } std::ifstream file(dictfile.c_str()); if (file.good()) { @@ -206,7 +205,7 @@ void run_trac_ca_its(bool cosmics = false, memParams[0].CellsMemoryCoefficients[iLayer] = 0.001f; } } else { - // Comment for pp + // PbPb tracking params // ---- trackParams.resize(3); memParams.resize(3); @@ -249,11 +248,7 @@ void run_trac_ca_its(bool cosmics = false, pattIt = patt.begin(); int rofId{0}; for (auto& rof : *rofs) { - counterEvent++; - // if (counterEvent != 111) { - // continue; - // } - std::cout << "ROF: " << counterEvent << std::endl; + auto start = std::chrono::high_resolution_clock::now(); auto it = pattIt; o2::its::ioutils::loadROFrameData(rof, event, clSpan, pattIt, dict, labels); @@ -275,9 +270,7 @@ void run_trac_ca_its(bool cosmics = false, if (!vertITS.empty()) { // Using only the first vertex in the list std::cout << " - Reconstructed vertex: x = " << vertITS[0].getX() << " y = " << vertITS[0].getY() << " x = " << vertITS[0].getZ() << std::endl; - // LOG(FATAL) << event.getPrimaryVertex(0).x << " " << event.getPrimaryVertex(0).y << " " << event.getPrimaryVertex(0).z; - // LOG(FATAL) << " --- " << event.getPrimaryVerticesNum(); - // event.addPrimaryVertex(vertITS[0].getX(), vertITS[0].getY(), vertITS[0].getZ()); + event.addPrimaryVertex(vertITS[0].getX(), vertITS[0].getY(), vertITS[0].getZ()); } else { std::cout << " - Vertex not reconstructed, tracking skipped" << std::endl; } From a6bb2452ddbefd2dc7963405c617e308e2459f65 Mon Sep 17 00:00:00 2001 From: Matteo Concas Date: Tue, 14 Sep 2021 14:08:28 +0200 Subject: [PATCH 06/10] Update Smoother.h --- .../ITS/tracking/include/ITStracking/Smoother.h | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Smoother.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Smoother.h index 6d075cead0421..e83251e368bd2 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Smoother.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Smoother.h @@ -1,13 +1,14 @@ -// 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". +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. // -// See http://alice-o2.web.cern.ch/license for full licensing information. +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". // // 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 Smoother.h /// \brief Class to handle Kalman smoothing for ITS tracking. /// Its instance stores the state of the track to the level we want to smooth to avoid multiple re-propagations when testing different clusters. @@ -65,4 +66,4 @@ class Smoother #endif }; } // namespace its -} // namespace o2 \ No newline at end of file +} // namespace o2 From 2aa0b40300bbe2d7004544ee58c949302caedf6d Mon Sep 17 00:00:00 2001 From: Matteo Concas Date: Tue, 14 Sep 2021 14:08:59 +0200 Subject: [PATCH 07/10] Update Smoother.cxx --- Detectors/ITSMFT/ITS/tracking/src/Smoother.cxx | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/src/Smoother.cxx b/Detectors/ITSMFT/ITS/tracking/src/Smoother.cxx index 6743ef2e865c1..04fe093af8d0f 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Smoother.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Smoother.cxx @@ -1,12 +1,15 @@ -// 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". +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. // -// See http://alice-o2.web.cern.ch/license for full licensing information. +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". // // 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. +// +// \author matteo.concas@cern.ch #include "ITStracking/Smoother.h" @@ -228,4 +231,4 @@ bool Smoother::testCluster(const int clusterId, const ROframe& event) template class Smoother<7>; } // namespace its -} // namespace o2 \ No newline at end of file +} // namespace o2 From bec838d3118ac2cf005191f4c7e5d31703993b23 Mon Sep 17 00:00:00 2001 From: Matteo Concas Date: Thu, 16 Sep 2021 14:20:07 +0200 Subject: [PATCH 08/10] Rename macro --- .../{CheckTrackerFakeClusters.C => CheckFakeClustersTracks.C} | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) rename Detectors/ITSMFT/ITS/macros/test/{CheckTrackerFakeClusters.C => CheckFakeClustersTracks.C} (96%) diff --git a/Detectors/ITSMFT/ITS/macros/test/CheckTrackerFakeClusters.C b/Detectors/ITSMFT/ITS/macros/test/CheckFakeClustersTracks.C similarity index 96% rename from Detectors/ITSMFT/ITS/macros/test/CheckTrackerFakeClusters.C rename to Detectors/ITSMFT/ITS/macros/test/CheckFakeClustersTracks.C index 0a5b62c1212d6..75050ffa0a752 100644 --- a/Detectors/ITSMFT/ITS/macros/test/CheckTrackerFakeClusters.C +++ b/Detectors/ITSMFT/ITS/macros/test/CheckFakeClustersTracks.C @@ -13,7 +13,7 @@ #include "SimulationDataFormat/MCTruthContainer.h" #endif -void CheckMCTracks(const std::string tracfile = "o2trac_its.root") +void CheckFakeClustersTracks(const std::string tracfile = "o2trac_its.root") { TFile* file1 = TFile::Open(tracfile.data()); TTree* recTree = (TTree*)file1->Get("o2sim"); From 5349e7f5b500949bc8abdb18e48550f54be1aabb Mon Sep 17 00:00:00 2001 From: Matteo Concas Date: Thu, 16 Sep 2021 20:30:50 +0200 Subject: [PATCH 09/10] Fix mc flag in setting in WF --- .../ITS/include/DataFormatsITS/TrackITS.h | 4 +- .../ITSMFT/ITS/macros/test/CMakeLists.txt | 2 +- Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx | 3 - .../ITSMFT/ITS/workflow/src/TrackerSpec.cxx | 3 - macro/run_trac_ca_its.C | 57 ++++++++++--------- 5 files changed, 34 insertions(+), 35 deletions(-) diff --git a/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h b/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h index 6968aa8ade4c4..34f0f79ccfdad 100644 --- a/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h +++ b/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h @@ -95,7 +95,7 @@ class TrackITS : public o2::track::TrackParCov const o2::track::TrackParCov& getParamOut() const { return mParamOut; } void setPattern(uint32_t p) { mPattern = p; } - int getPattern() const { return mPattern; } + uint32_t getPattern() const { return mPattern; } bool hasHitOnLayer(int i) { return mPattern & (0x1 << i); } bool isFakeOnLayer(int i) { return !(mPattern & (0x1 << (16 + i))); } int getNFakeClusters(); @@ -133,6 +133,8 @@ class TrackITSExt : public TrackITS setNumberOfClusters(ncl); } + GPUdDefault() TrackITSExt(const TrackITSExt& t) = default; + void setClusterIndex(int l, int i) { int ncl = getNumberOfClusters(); diff --git a/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt b/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt index dcd6026bd9094..24f196e47531c 100644 --- a/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt @@ -70,7 +70,7 @@ o2_add_test_root_macro(RunGPUTracking.C PUBLIC_LINK_LIBRARIES O2::GPUTracking LABELS its) -o2_add_test_root_macro(CheckTrackerFakeClusters.C +o2_add_test_root_macro(CheckFakeClustersTracks.C PUBLIC_LINK_LIBRARIES O2::SimulationDataFormat O2::ITSBase O2::DataFormatsITS diff --git a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx index d5dbe3a2c353a..9f86c6f487af8 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx @@ -527,9 +527,6 @@ void Tracker::computeRoadsMClabels() void Tracker::computeTracksMClabels() { - if (!mTimeFrame->hasMCinformation()) { - return; - } for (int iROF{0}; iROF < mTimeFrame->getNrof(); ++iROF) { for (auto& track : mTimeFrame->getTracks(iROF)) { diff --git a/Detectors/ITSMFT/ITS/workflow/src/TrackerSpec.cxx b/Detectors/ITSMFT/ITS/workflow/src/TrackerSpec.cxx index b8b6cc9350bd0..c42c52f1cdea8 100644 --- a/Detectors/ITSMFT/ITS/workflow/src/TrackerSpec.cxx +++ b/Detectors/ITSMFT/ITS/workflow/src/TrackerSpec.cxx @@ -252,17 +252,14 @@ void TrackerDPL::run(ProcessingContext& pc) auto& trc{tracks[iTrk]}; trc.setFirstClusterEntry(allClusIdx.size()); // before adding tracks, create final cluster indices int ncl = trc.getNumberOfClusters(), nclf = 0; - uint8_t patt = 0; for (int ic = TrackITSExt::MaxClusters; ic--;) { // track internally keeps in->out cluster indices, but we want to store the references as out->in!!! auto clid = trc.getClusterIndex(ic); if (clid >= 0) { allClusIdx.push_back(clid); nclf++; - patt |= 0x1 << ic; } } assert(ncl == nclf); - trc.setPattern(patt); allTracks.emplace_back(trc); } } diff --git a/macro/run_trac_ca_its.C b/macro/run_trac_ca_its.C index 129adba7bfaf5..e4bcc887ad0ba 100644 --- a/macro/run_trac_ca_its.C +++ b/macro/run_trac_ca_its.C @@ -208,35 +208,37 @@ void run_trac_ca_its(bool cosmics = false, // PbPb tracking params // ---- trackParams.resize(3); - memParams.resize(3); trackParams[0].TrackletMaxDeltaPhi = 0.05f; - trackParams[1].TrackletMaxDeltaPhi = 0.1f; - trackParams[2].MinTrackLength = 5; - trackParams[2].TrackletMaxDeltaPhi = 0.3; + trackParams[0].DeltaROF = 0; + trackParams[1].CopyCuts(trackParams[0], 2.); + trackParams[1].DeltaROF = 0; + trackParams[2].CopyCuts(trackParams[1], 2.); + trackParams[2].DeltaROF = 1; + trackParams[2].MinTrackLength = 4; + memParams.resize(3); // --- // Uncomment for pp - trackParams.resize(2); - std::array kmaxDCAxy1 = {1.f * 2.0, 0.4f * 2.0, 0.4f * 2.0, 2.0f * 2.0, 3.f * 2.0}; - std::array kmaxDCAz1 = {1.f * 2.0, 0.4f * 2.0, 0.4f * 2.0, 2.0f * 2.0, 3.f * 2.0}; - std::array kmaxDN1 = {0.005f * 2.0, 0.0035f * 2.0, 0.009f * 2.0, 0.03f * 2.0}; - std::array kmaxDP1 = {0.02f * 2.0, 0.005f * 2.0, 0.006f * 2.0, 0.007f * 2.0}; - std::array kmaxDZ1 = {1.f * 2.0, 1.f * 2.0, 2.0f * 2.0, 2.0f * 2.0, 2.0f * 2.0, 2.0f * 2.0}; - const float kDoublTanL1 = 0.05f * 5.; - const float kDoublPhi1 = 0.2f * 5.; - trackParams[1].MinTrackLength = 4; - trackParams[1].TrackletMaxDeltaPhi = 0.3; - trackParams[1].CellMaxDeltaPhi = 0.2 * 2; - trackParams[1].CellMaxDeltaTanLambda = 0.05 * 2; - std::copy(kmaxDZ1.begin(), kmaxDZ1.end(), trackParams[1].TrackletMaxDeltaZ.begin()); - std::copy(kmaxDCAxy1.begin(), kmaxDCAxy1.end(), trackParams[1].CellMaxDCA.begin()); - std::copy(kmaxDCAz1.begin(), kmaxDCAz1.end(), trackParams[1].CellMaxDeltaZ.begin()); - std::copy(kmaxDP1.begin(), kmaxDP1.end(), trackParams[1].NeighbourMaxDeltaCurvature.begin()); - std::copy(kmaxDN1.begin(), kmaxDN1.end(), trackParams[1].NeighbourMaxDeltaN.begin()); - memParams.resize(2); + // trackParams.resize(2); + // std::array kmaxDCAxy1 = {1.f * 2.0, 0.4f * 2.0, 0.4f * 2.0, 2.0f * 2.0, 3.f * 2.0}; + // std::array kmaxDCAz1 = {1.f * 2.0, 0.4f * 2.0, 0.4f * 2.0, 2.0f * 2.0, 3.f * 2.0}; + // std::array kmaxDN1 = {0.005f * 2.0, 0.0035f * 2.0, 0.009f * 2.0, 0.03f * 2.0}; + // std::array kmaxDP1 = {0.02f * 2.0, 0.005f * 2.0, 0.006f * 2.0, 0.007f * 2.0}; + // std::array kmaxDZ1 = {1.f * 2.0, 1.f * 2.0, 2.0f * 2.0, 2.0f * 2.0, 2.0f * 2.0, 2.0f * 2.0}; + // const float kDoublTanL1 = 0.05f * 5.; + // const float kDoublPhi1 = 0.2f * 5.; + // trackParams[1].MinTrackLength = 4; + // trackParams[1].TrackletMaxDeltaPhi = 0.3; + // trackParams[1].CellMaxDeltaPhi = 0.2 * 2; + // trackParams[1].CellMaxDeltaTanLambda = 0.05 * 2; + // std::copy(kmaxDZ1.begin(), kmaxDZ1.end(), trackParams[1].TrackletMaxDeltaZ.begin()); + // std::copy(kmaxDCAxy1.begin(), kmaxDCAxy1.end(), trackParams[1].CellMaxDCA.begin()); + // std::copy(kmaxDCAz1.begin(), kmaxDCAz1.end(), trackParams[1].CellMaxDeltaZ.begin()); + // std::copy(kmaxDP1.begin(), kmaxDP1.end(), trackParams[1].NeighbourMaxDeltaCurvature.begin()); + // std::copy(kmaxDN1.begin(), kmaxDN1.end(), trackParams[1].NeighbourMaxDeltaN.begin()); + // memParams.resize(2); // --- } - int currentEvent = -1; gsl::span patt(patterns->data(), patterns->size()); auto pattIt = patt.begin(); @@ -286,10 +288,6 @@ void run_trac_ca_its(bool cosmics = false, tf.printVertices(); o2::its::Tracker tracker(new o2::its::TrackerTraitsCPU(&tf)); - tracker.setBz(field->getBz(origD)); - tracker.setParameters(memParams, trackParams); - tracker.clustersToTracks(); - //-------- init lookuptable --------// if (useLUT) { auto* lut = o2::base::MatLayerCylSet::loadFromFile(matLUTFile); o2::base::Propagator::Instance()->setMatLUT(lut); @@ -303,6 +301,11 @@ void run_trac_ca_its(bool cosmics = false, LOG(INFO) << "Material LUT " << matLUTFile << " file is absent, only TGeo can be used"; } + tracker.setBz(field->getBz(origD)); + tracker.setParameters(memParams, trackParams); + tracker.clustersToTracks(); + //-------- init lookuptable --------// + for (int iROF{0}; iROF < tf.getNrof(); ++iROF) { tracksITS.clear(); for (auto& trc : tf.getTracks(iROF)) { From 9f51b428c4d304c648e6c78e90a2aaeb53b27cd0 Mon Sep 17 00:00:00 2001 From: Matteo Concas Date: Fri, 17 Sep 2021 14:42:09 +0200 Subject: [PATCH 10/10] Extend ChecTracksCA.C --- .../ITSMFT/ITS/macros/test/CMakeLists.txt | 7 -- .../ITS/macros/test/CheckFakeClustersTracks.C | 69 -------------- .../ITSMFT/ITS/macros/test/CheckTracksCA.C | 94 ++++++++++++++++++- macro/run_trac_ca_its.C | 54 +++++------ 4 files changed, 120 insertions(+), 104 deletions(-) delete mode 100644 Detectors/ITSMFT/ITS/macros/test/CheckFakeClustersTracks.C diff --git a/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt b/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt index 24f196e47531c..befc776289e25 100644 --- a/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt @@ -70,13 +70,6 @@ o2_add_test_root_macro(RunGPUTracking.C PUBLIC_LINK_LIBRARIES O2::GPUTracking LABELS its) -o2_add_test_root_macro(CheckFakeClustersTracks.C - PUBLIC_LINK_LIBRARIES O2::SimulationDataFormat - O2::ITSBase - O2::DataFormatsITS - O2::DataFormatsITSMFT - LABELS its) - o2_add_test_root_macro(DisplayTrack.C PUBLIC_LINK_LIBRARIES O2::ITSBase O2::DataFormatsITSMFT diff --git a/Detectors/ITSMFT/ITS/macros/test/CheckFakeClustersTracks.C b/Detectors/ITSMFT/ITS/macros/test/CheckFakeClustersTracks.C deleted file mode 100644 index 75050ffa0a752..0000000000000 --- a/Detectors/ITSMFT/ITS/macros/test/CheckFakeClustersTracks.C +++ /dev/null @@ -1,69 +0,0 @@ -#if !defined(__CLING__) || defined(__ROOTCLING__) -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include "DataFormatsITS/TrackITS.h" -#include "SimulationDataFormat/MCCompLabel.h" -#include "SimulationDataFormat/MCTruthContainer.h" -#endif - -void CheckFakeClustersTracks(const std::string tracfile = "o2trac_its.root") -{ - TFile* file1 = TFile::Open(tracfile.data()); - TTree* recTree = (TTree*)file1->Get("o2sim"); - std::vector* recArr = nullptr; - recTree->SetBranchAddress("ITSTrack", &recArr); - std::vector* trkLabArr = nullptr; - recTree->SetBranchAddress("ITSTrackMCTruth", &trkLabArr); - std::vector histLength; - std::vector histLength1Fake; - histLength.resize(4); - histLength1Fake.resize(4); - for (int iH{4}; iH < 8; ++iH) { - histLength[iH - 4] = new TH1I(Form("trk_len_%d", iH), Form("trk_len = %d", iH), 7, -.5, 6.5); - } - for (int iH{4}; iH < 8; ++iH) { - histLength1Fake[iH - 4] = new TH1I(Form("trk_len_%d_1f", iH), Form("trk_len = %d, 1 Fake", iH), 7, -.5, 6.5); - } - auto nt = new TNtuple("TrackITSInfo", "Info about ITS tracks", "len:l0:l1:l2:l3:l4:l5:6"); - std::array labels; - for (unsigned int iEntry{0}; iEntry < recTree->GetEntriesFast(); ++iEntry) { - recTree->GetEntry(iEntry); - for (unsigned int iTrack{0}; iTrack < recArr->size(); ++iTrack) { - auto& track = (*recArr)[iTrack]; - auto len = (*recArr)[iTrack].getNClusters(); - int count = 0; - for (int iLabel{0}; iLabel < 7; ++iLabel) { - if (track.hasHitOnLayer(iLabel)) { - count++; - if (track.isFakeOnLayer(iLabel)) { - histLength[len - 4]->Fill(iLabel); - if (track.getNFakeClusters() == 1) { - histLength1Fake[len - 4]->Fill(iLabel); - } - } - } - } - if (count != len) { - std::cout << "mismatch!\n"; - } - } - } - auto canvas = new TCanvas("c1", "c1", 2400, 1600); - canvas->Divide(4, 2); - for (int iH{0}; iH < 4; ++iH) { - canvas->cd(iH + 1); - histLength[iH]->Draw(); - } - for (int iH{0}; iH < 4; ++iH) { - canvas->cd(iH + 5); - histLength1Fake[iH]->Draw(); - } - canvas->SaveAs("fakeClusters.png", "recreate"); -} \ No newline at end of file diff --git a/Detectors/ITSMFT/ITS/macros/test/CheckTracksCA.C b/Detectors/ITSMFT/ITS/macros/test/CheckTracksCA.C index ef37ede67321c..e6b4b91e00c6f 100644 --- a/Detectors/ITSMFT/ITS/macros/test/CheckTracksCA.C +++ b/Detectors/ITSMFT/ITS/macros/test/CheckTracksCA.C @@ -12,6 +12,9 @@ #include #include #include +#include +#include +#include #include "ITSBase/GeometryTGeo.h" #include "SimulationDataFormat/TrackReference.h" @@ -43,7 +46,7 @@ struct ParticleInfo { #pragma link C++ class ParticleInfo + ; -void CheckTracksCA(std::string tracfile = "o2trac_its.root", std::string clusfile = "o2clus_its.root", std::string kinefile = "o2sim_Kine.root") +void CheckTracksCA(bool doFakeClStud = true, std::string tracfile = "o2trac_its.root", std::string clusfile = "o2clus_its.root", std::string kinefile = "o2sim_Kine.root") { using namespace o2::itsmft; @@ -271,4 +274,93 @@ void CheckTracksCA(std::string tracfile = "o2trac_its.root", std::string clusfil clone->Write("clones"); file.Close(); std::cout << " done." << std::endl; + + ////////////////////// + // Fake clusters study + if (doFakeClStud) { + std::vector histLength, histLength1Fake, histLengthNoCl, histLength1FakeNoCl; + std::vector stackLength, stackLength1Fake; + std::vector legends, legends1Fake; + histLength.resize(4); + histLength1Fake.resize(4); + histLengthNoCl.resize(4); + histLength1FakeNoCl.resize(4); + stackLength.resize(4); + stackLength1Fake.resize(4); + legends.resize(4); + legends1Fake.resize(4); + + for (int iH{4}; iH < 8; ++iH) { + histLength[iH - 4] = new TH1I(Form("trk_len_%d", iH), "#exists cluster", 7, -.5, 6.5); + histLength[iH - 4]->SetFillColor(kBlue); + histLength[iH - 4]->SetLineColor(kBlue); + histLength[iH - 4]->SetFillStyle(3352); + histLengthNoCl[iH - 4] = new TH1I(Form("trk_len_%d_nocl", iH), "#slash{#exists} cluster", 7, -.5, 6.5); + histLengthNoCl[iH - 4]->SetFillColor(kRed); + histLengthNoCl[iH - 4]->SetLineColor(kRed); + histLengthNoCl[iH - 4]->SetFillStyle(3352); + stackLength[iH - 4] = new THStack(Form("stack_trk_len_%d", iH), Form("trk_len=%d", iH)); + stackLength[iH - 4]->Add(histLength[iH - 4]); + stackLength[iH - 4]->Add(histLengthNoCl[iH - 4]); + } + for (int iH{4}; iH < 8; ++iH) { + histLength1Fake[iH - 4] = new TH1I(Form("trk_len_%d_1f", iH), "#exists cluster", 7, -.5, 6.5); + histLength1Fake[iH - 4]->SetFillColor(kBlue); + histLength1Fake[iH - 4]->SetLineColor(kBlue); + histLength1Fake[iH - 4]->SetFillStyle(3352); + histLength1FakeNoCl[iH - 4] = new TH1I(Form("trk_len_%d_1f_nocl", iH), "#slash{#exists} cluster", 7, -.5, 6.5); + histLength1FakeNoCl[iH - 4]->SetFillColor(kRed); + histLength1FakeNoCl[iH - 4]->SetLineColor(kRed); + histLength1FakeNoCl[iH - 4]->SetFillStyle(3352); + stackLength1Fake[iH - 4] = new THStack(Form("stack_trk_len_%d_1f", iH), Form("trk_len=%d, 1 Fake", iH)); + stackLength1Fake[iH - 4]->Add(histLength1Fake[iH - 4]); + stackLength1Fake[iH - 4]->Add(histLength1FakeNoCl[iH - 4]); + } + + for (auto& event : info) { + for (auto& part : event) { + int nCl{0}; + for (unsigned int bit{0}; bit < sizeof(pInfo.clusters) * 8; ++bit) { + nCl += bool(part.clusters & (1 << bit)); + } + if (nCl < 3) { + continue; + } + + auto& track = part.track; + auto len = track.getNClusters(); + for (int iLayer{0}; iLayer < 7; ++iLayer) { + if (track.hasHitOnLayer(iLayer)) { + if (track.isFakeOnLayer(iLayer)) { // Reco track has fake cluster + if (part.clusters & (0x1 << iLayer)) { // Correct cluster exists + histLength[len - 4]->Fill(iLayer); + if (track.getNFakeClusters() == 1) { + histLength1Fake[len - 4]->Fill(iLayer); + } + } else { + histLengthNoCl[len - 4]->Fill(iLayer); + if (track.getNFakeClusters() == 1) { + histLength1FakeNoCl[len - 4]->Fill(iLayer); + } + } + } + } + } + } + } + + auto canvas = new TCanvas("fc_canvas", "Fake clusters", 1600, 1000); + canvas->Divide(4, 2); + for (int iH{0}; iH < 4; ++iH) { + canvas->cd(iH + 1); + stackLength[iH]->Draw(); + gPad->BuildLegend(); + } + for (int iH{0}; iH < 4; ++iH) { + canvas->cd(iH + 5); + stackLength1Fake[iH]->Draw(); + gPad->BuildLegend(); + } + canvas->SaveAs("fakeClusters.png", "recreate"); + } } diff --git a/macro/run_trac_ca_its.C b/macro/run_trac_ca_its.C index e4bcc887ad0ba..90c99753f8e18 100644 --- a/macro/run_trac_ca_its.C +++ b/macro/run_trac_ca_its.C @@ -207,35 +207,35 @@ void run_trac_ca_its(bool cosmics = false, } else { // PbPb tracking params // ---- - trackParams.resize(3); - trackParams[0].TrackletMaxDeltaPhi = 0.05f; - trackParams[0].DeltaROF = 0; - trackParams[1].CopyCuts(trackParams[0], 2.); - trackParams[1].DeltaROF = 0; - trackParams[2].CopyCuts(trackParams[1], 2.); - trackParams[2].DeltaROF = 1; - trackParams[2].MinTrackLength = 4; - memParams.resize(3); + // trackParams.resize(3); + // trackParams[0].TrackletMaxDeltaPhi = 0.05f; + // trackParams[0].DeltaROF = 0; + // trackParams[1].CopyCuts(trackParams[0], 2.); + // trackParams[1].DeltaROF = 0; + // trackParams[2].CopyCuts(trackParams[1], 2.); + // trackParams[2].DeltaROF = 1; + // trackParams[2].MinTrackLength = 4; + // memParams.resize(3); // --- // Uncomment for pp - // trackParams.resize(2); - // std::array kmaxDCAxy1 = {1.f * 2.0, 0.4f * 2.0, 0.4f * 2.0, 2.0f * 2.0, 3.f * 2.0}; - // std::array kmaxDCAz1 = {1.f * 2.0, 0.4f * 2.0, 0.4f * 2.0, 2.0f * 2.0, 3.f * 2.0}; - // std::array kmaxDN1 = {0.005f * 2.0, 0.0035f * 2.0, 0.009f * 2.0, 0.03f * 2.0}; - // std::array kmaxDP1 = {0.02f * 2.0, 0.005f * 2.0, 0.006f * 2.0, 0.007f * 2.0}; - // std::array kmaxDZ1 = {1.f * 2.0, 1.f * 2.0, 2.0f * 2.0, 2.0f * 2.0, 2.0f * 2.0, 2.0f * 2.0}; - // const float kDoublTanL1 = 0.05f * 5.; - // const float kDoublPhi1 = 0.2f * 5.; - // trackParams[1].MinTrackLength = 4; - // trackParams[1].TrackletMaxDeltaPhi = 0.3; - // trackParams[1].CellMaxDeltaPhi = 0.2 * 2; - // trackParams[1].CellMaxDeltaTanLambda = 0.05 * 2; - // std::copy(kmaxDZ1.begin(), kmaxDZ1.end(), trackParams[1].TrackletMaxDeltaZ.begin()); - // std::copy(kmaxDCAxy1.begin(), kmaxDCAxy1.end(), trackParams[1].CellMaxDCA.begin()); - // std::copy(kmaxDCAz1.begin(), kmaxDCAz1.end(), trackParams[1].CellMaxDeltaZ.begin()); - // std::copy(kmaxDP1.begin(), kmaxDP1.end(), trackParams[1].NeighbourMaxDeltaCurvature.begin()); - // std::copy(kmaxDN1.begin(), kmaxDN1.end(), trackParams[1].NeighbourMaxDeltaN.begin()); - // memParams.resize(2); + trackParams.resize(2); + std::array kmaxDCAxy1 = {1.f * 2.0, 0.4f * 2.0, 0.4f * 2.0, 2.0f * 2.0, 3.f * 2.0}; + std::array kmaxDCAz1 = {1.f * 2.0, 0.4f * 2.0, 0.4f * 2.0, 2.0f * 2.0, 3.f * 2.0}; + std::array kmaxDN1 = {0.005f * 2.0, 0.0035f * 2.0, 0.009f * 2.0, 0.03f * 2.0}; + std::array kmaxDP1 = {0.02f * 2.0, 0.005f * 2.0, 0.006f * 2.0, 0.007f * 2.0}; + std::array kmaxDZ1 = {1.f * 2.0, 1.f * 2.0, 2.0f * 2.0, 2.0f * 2.0, 2.0f * 2.0, 2.0f * 2.0}; + const float kDoublTanL1 = 0.05f * 5.; + const float kDoublPhi1 = 0.2f * 5.; + trackParams[1].MinTrackLength = 4; + trackParams[1].TrackletMaxDeltaPhi = 0.3; + trackParams[1].CellMaxDeltaPhi = 0.2 * 2; + trackParams[1].CellMaxDeltaTanLambda = 0.05 * 2; + std::copy(kmaxDZ1.begin(), kmaxDZ1.end(), trackParams[1].TrackletMaxDeltaZ.begin()); + std::copy(kmaxDCAxy1.begin(), kmaxDCAxy1.end(), trackParams[1].CellMaxDCA.begin()); + std::copy(kmaxDCAz1.begin(), kmaxDCAz1.end(), trackParams[1].CellMaxDeltaZ.begin()); + std::copy(kmaxDP1.begin(), kmaxDP1.end(), trackParams[1].NeighbourMaxDeltaCurvature.begin()); + std::copy(kmaxDN1.begin(), kmaxDN1.end(), trackParams[1].NeighbourMaxDeltaN.begin()); + memParams.resize(2); // --- }