diff --git a/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h b/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h index d4eb494d4dfa4..62bf846f65205 100644 --- a/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h +++ b/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h @@ -18,6 +18,7 @@ #include +#include "GPUCommonDef.h" #include "ReconstructionDataFormats/Track.h" #include "CommonDataFormat/RangeReference.h" @@ -33,6 +34,10 @@ namespace its class TrackITS : public o2::track::TrackParCov { + enum UserBits { + kNextROF = 1 + }; + using Cluster = o2::itsmft::Cluster; using ClusRefs = o2::dataformats::RangeRefComp<4>; @@ -40,13 +45,13 @@ class TrackITS : public o2::track::TrackParCov using o2::track::TrackParCov::TrackParCov; // inherit base constructors static constexpr int MaxClusters = 16; - TrackITS() = default; - TrackITS(const TrackITS& t) = default; - TrackITS(const o2::track::TrackParCov& parcov) : o2::track::TrackParCov{parcov} {} - TrackITS(const o2::track::TrackParCov& parCov, float chi2, const o2::track::TrackParCov& outer) + GPUdDefault() TrackITS() = default; + GPUdDefault() TrackITS(const TrackITS& t) = default; + GPUd() TrackITS(const o2::track::TrackParCov& parcov) : o2::track::TrackParCov{parcov} {} + GPUd() TrackITS(const o2::track::TrackParCov& parCov, float chi2, const o2::track::TrackParCov& outer) : o2::track::TrackParCov{parCov}, mParamOut{outer}, mChi2{chi2} {} - TrackITS& operator=(const TrackITS& tr) = default; - ~TrackITS() = default; + GPUdDefault() TrackITS& operator=(const TrackITS& tr) = default; + GPUdDefault() ~TrackITS() = default; // These functions must be provided bool propagate(float alpha, float x, float bz); @@ -93,6 +98,9 @@ class TrackITS : public o2::track::TrackParCov int getPattern() const { return mPattern; } bool hasHitOnLayer(int i) { return mPattern & (0x1 << i); } + void setNextROFbit(bool toggle = true) { setUserField((getUserField() & ~kNextROF) | (-toggle & kNextROF)); } + bool hasHitInNextROF() const { return getUserField() & kNextROF; } + private: o2::track::TrackParCov mParamOut; ///< parameter at largest radius ClusRefs mClusRef; ///< references on clusters @@ -109,8 +117,8 @@ class TrackITSExt : public TrackITS static constexpr int MaxClusters = 16; /// Prepare for overlaps and new detector configurations using TrackITS::TrackITS; // inherit base constructors - TrackITSExt(o2::track::TrackParCov&& parCov, short ncl, float chi2, - o2::track::TrackParCov&& outer, std::array cls) + GPUd() TrackITSExt(o2::track::TrackParCov&& parCov, short ncl, float chi2, + o2::track::TrackParCov&& outer, std::array cls) : TrackITS(parCov, chi2, outer), mIndex{cls} { setNumberOfClusters(ncl); diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h index 7ab2c6d304686..c23cb7e7d1570 100644 --- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h +++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h @@ -213,8 +213,8 @@ class TrackParametrization GPUd() bool isValid() const; GPUd() void invalidate(); - GPUd() uint16_t getUserField() const; - GPUd() void setUserField(uint16_t v); + GPUhd() uint16_t getUserField() const; + GPUhd() void setUserField(uint16_t v); GPUd() void printParam() const; #ifndef GPUCA_ALIGPUCODE @@ -655,13 +655,13 @@ GPUdi() void TrackParametrization::invalidate() } template -GPUdi() uint16_t TrackParametrization::getUserField() const +GPUhdi() uint16_t TrackParametrization::getUserField() const { return mUserField; } template -GPUdi() void TrackParametrization::setUserField(uint16_t v) +GPUhdi() void TrackParametrization::setUserField(uint16_t v) { mUserField = v; } diff --git a/DataFormats/simulation/include/SimulationDataFormat/MCCompLabel.h b/DataFormats/simulation/include/SimulationDataFormat/MCCompLabel.h index 6436958e7bb34..ec81285c9e802 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/MCCompLabel.h +++ b/DataFormats/simulation/include/SimulationDataFormat/MCCompLabel.h @@ -21,13 +21,13 @@ namespace o2 class MCCompLabel { private: - static constexpr ULong64_t ul0x1 = 0x1; - static constexpr ULong64_t NotSet = 0xffffffffffffffff; - static constexpr ULong64_t Noise = 0xfffffffffffffffe; - static constexpr ULong64_t Fake = ul0x1 << 63; + static constexpr uint64_t ul0x1 = 0x1; + static constexpr uint64_t NotSet = 0xffffffffffffffff; + static constexpr uint64_t Noise = 0xfffffffffffffffe; + static constexpr uint64_t Fake = ul0x1 << 63; static constexpr int NReservedBits = 1; - ULong64_t mLabel = NotSet; ///< MC label encoding MCtrack ID and MCevent origin + uint64_t mLabel = NotSet; ///< MC label encoding MCtrack ID and MCevent origin public: // number of bits reserved for MC track ID, DON'T modify this, since the @@ -38,17 +38,17 @@ class MCCompLabel // the rest of the bits is reserved at the moment // check if the fields are defined consistently - static_assert(nbitsTrackID + nbitsEvID + nbitsSrcID <= sizeof(ULong64_t) * 8 - NReservedBits, + static_assert(nbitsTrackID + nbitsEvID + nbitsSrcID <= sizeof(uint64_t) * 8 - NReservedBits, "Fields cannot be stored in 64 bits"); // mask to extract MC track ID - static constexpr ULong64_t maskTrackID = (ul0x1 << nbitsTrackID) - 1; + static constexpr uint64_t maskTrackID = (ul0x1 << nbitsTrackID) - 1; // mask to extract MC track ID - static constexpr ULong64_t maskEvID = (ul0x1 << nbitsEvID) - 1; + static constexpr uint64_t maskEvID = (ul0x1 << nbitsEvID) - 1; // mask to extract MC track ID - static constexpr ULong64_t maskSrcID = (ul0x1 << nbitsSrcID) - 1; + static constexpr uint64_t maskSrcID = (ul0x1 << nbitsSrcID) - 1; // mask for all used fields - static constexpr ULong64_t maskFull = (ul0x1 << (nbitsTrackID + nbitsEvID + nbitsSrcID)) - 1; + static constexpr uint64_t maskFull = (ul0x1 << (nbitsTrackID + nbitsEvID + nbitsSrcID)) - 1; MCCompLabel(int trackID, int evID, int srcID, bool fake = false) { set(trackID, evID, srcID, fake); } MCCompLabel(bool noise = false) @@ -88,7 +88,7 @@ class MCCompLabel } // allow to retrieve bare label - ULong64_t getRawValue() const { return mLabel; } + uint64_t getRawValue() const { return mLabel; } // comparison operator, compares only label, not eventual weight or correctness info bool operator==(const MCCompLabel& other) const { return (mLabel & maskFull) == (other.mLabel & maskFull); } @@ -112,9 +112,9 @@ class MCCompLabel void set(unsigned int trackID, int evID, int srcID, bool fake) { /// compose label: the track 1st cast to UInt32_t to preserve the sign! - mLabel = (maskTrackID & static_cast(trackID)) | - (maskEvID & static_cast(evID)) << nbitsTrackID | - (maskSrcID & static_cast(srcID)) << (nbitsTrackID + nbitsEvID); + mLabel = (maskTrackID & static_cast(trackID)) | + (maskEvID & static_cast(evID)) << nbitsTrackID | + (maskSrcID & static_cast(srcID)) << (nbitsTrackID + nbitsEvID); if (fake) { setFakeFlag(); } @@ -124,7 +124,7 @@ class MCCompLabel int getTrackIDSigned() const { return isFake() ? -getTrackID() : getTrackID(); } int getEventID() const { return (mLabel >> nbitsTrackID) & maskEvID; } int getSourceID() const { return (mLabel >> (nbitsTrackID + nbitsEvID)) & maskSrcID; } - ULong64_t getTrackEventSourceID() const { return static_cast(mLabel & maskFull); } + uint64_t getTrackEventSourceID() const { return static_cast(mLabel & maskFull); } void get(int& trackID, int& evID, int& srcID, bool& fake) { /// parse label diff --git a/DataFormats/simulation/include/SimulationDataFormat/MCTruthContainer.h b/DataFormats/simulation/include/SimulationDataFormat/MCTruthContainer.h index 83cd59e58db97..92156bd35db78 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/MCTruthContainer.h +++ b/DataFormats/simulation/include/SimulationDataFormat/MCTruthContainer.h @@ -43,7 +43,7 @@ struct MCTruthHeaderElement { MCTruthHeaderElement() = default; // for ROOT IO MCTruthHeaderElement(uint32_t i) : index(i) {} - uint32_t index = -1; // the index into the actual MC track storage (-1 if invalid) + uint32_t index = (uint32_t)-1; // the index into the actual MC track storage (-1 if invalid) ClassDefNV(MCTruthHeaderElement, 1); }; diff --git a/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt b/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt index 064ebcbef14b4..6894e9fc8b147 100644 --- a/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt @@ -59,6 +59,13 @@ o2_add_test_root_macro(CheckTracks.C O2::DataFormatsITSMFT LABELS its) +o2_add_test_root_macro(CheckTracksCA.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/CheckTracks.C b/Detectors/ITSMFT/ITS/macros/test/CheckTracks.C index 27f382b289e36..a9374991592fc 100644 --- a/Detectors/ITSMFT/ITS/macros/test/CheckTracks.C +++ b/Detectors/ITSMFT/ITS/macros/test/CheckTracks.C @@ -153,14 +153,16 @@ void CheckTracks(std::string tracfile = "o2trac_its.root", std::string clusfile if (!recTree->GetEvent(frame)) continue; int loadedEventTracks = frame; + cout << "Number of tracks in frame " << frame << ": " << recArr->size() << std::endl; for (unsigned int i = 0; i < recArr->size(); i++) { // Find the last MC event within this reconstructed entry auto lab = (*trkLabArr)[i]; if (!lab.isValid()) { const TrackITS& recTrack = (*recArr)[i]; fak->Fill(recTrack.getPt()); } - if (!lab.isValid() || lab.getSourceID() != 0 || lab.getEventID() < 0 || lab.getEventID() >= nev) + if (!lab.isValid() || lab.getSourceID() != 0 || lab.getEventID() < 0 || lab.getEventID() >= nev) { continue; + } trackFrames[lab.getEventID()].update(frame, i); } } diff --git a/Detectors/ITSMFT/ITS/macros/test/CheckTracksCA.C b/Detectors/ITSMFT/ITS/macros/test/CheckTracksCA.C new file mode 100644 index 0000000000000..ef37ede67321c --- /dev/null +++ b/Detectors/ITSMFT/ITS/macros/test/CheckTracksCA.C @@ -0,0 +1,274 @@ +/// \file CheckTracksCA.C +/// \brief Simple macro to check ITSU tracks + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "ITSBase/GeometryTGeo.h" +#include "SimulationDataFormat/TrackReference.h" +#include "SimulationDataFormat/MCTrack.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsITS/TrackITS.h" + +#endif + +using namespace std; + +struct ParticleInfo { + int event; + int pdg; + float pt; + float eta; + float phi; + int mother; + int first; + unsigned short clusters = 0u; + unsigned char isReco = 0u; + unsigned char isFake = 0u; + bool isPrimary = 0u; + unsigned char storedStatus = 2; /// not stored = 2, fake = 1, good = 0 + o2::its::TrackITS track; +}; + +#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") +{ + + using namespace o2::itsmft; + using namespace o2::its; + + // Geometry + o2::base::GeometryManager::loadGeometry(); + auto gman = o2::its::GeometryTGeo::Instance(); + + // MC tracks + TFile* file0 = TFile::Open(kinefile.data()); + TTree* mcTree = (TTree*)gFile->Get("o2sim"); + mcTree->SetBranchStatus("*", 0); //disable all branches + mcTree->SetBranchStatus("MCTrack*", 1); + + std::vector* mcArr = nullptr; + mcTree->SetBranchAddress("MCTrack", &mcArr); + + // Clusters + TFile::Open(clusfile.data()); + TTree* clusTree = (TTree*)gFile->Get("o2sim"); + std::vector* clusArr = nullptr; + clusTree->SetBranchAddress("ITSClusterComp", &clusArr); + + // Cluster MC labels + o2::dataformats::MCTruthContainer* clusLabArr = nullptr; + clusTree->SetBranchAddress("ITSClusterMCTruth", &clusLabArr); + + // 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); + + std::cout << "** Filling particle table ... " << std::flush; + int lastEventIDcl = -1, cf = 0; + int nev = mcTree->GetEntriesFast(); + std::vector> info(nev); + for (int n = 0; n < nev; n++) { // loop over MC events + mcTree->GetEvent(n); + info[n].resize(mcArr->size()); + for (unsigned int mcI{0}; mcI < mcArr->size(); ++mcI) { + auto part = mcArr->at(mcI); + info[n][mcI].event = n; + info[n][mcI].pdg = part.GetPdgCode(); + info[n][mcI].pt = part.GetPt(); + info[n][mcI].phi = part.GetPhi(); + info[n][mcI].eta = part.GetEta(); + info[n][mcI].isPrimary = part.isPrimary(); + } + } + std::cout << "done." << std::endl; + + std::cout << "** Creating particle/clusters correspondance ... " << std::flush; + for (int frame = 0; frame < clusTree->GetEntriesFast(); frame++) { // Cluster frames + if (!clusTree->GetEvent(frame)) + continue; + + for (unsigned int iClus{0}; iClus < clusArr->size(); ++iClus) { + auto lab = (clusLabArr->getLabels(iClus))[0]; + if (!lab.isValid() || lab.getSourceID() != 0 || !lab.isCorrect()) + continue; + + int trackID, evID, srcID; + bool fake; + lab.get(trackID, evID, srcID, fake); + if (evID < 0 || evID >= (int)info.size()) { + std::cout << "Cluster MC label eventID out of range" << std::endl; + continue; + } + if (trackID < 0 || trackID >= (int)info[evID].size()) { + std::cout << "Cluster MC label trackID out of range" << std::endl; + continue; + } + + const CompClusterExt& c = (*clusArr)[iClus]; + auto layer = gman->getLayer(c.getSensorID()); + info[evID][trackID].clusters |= 1 << layer; + } + } + std::cout << "done." << std::endl; + + std::cout << "** Analysing tracks ... " << std::flush; + int unaccounted{0}, good{0}, fakes{0}, total{0}; + for (int frame = 0; frame < recTree->GetEntriesFast(); frame++) { // Cluster frames + if (!recTree->GetEvent(frame)) + continue; + total += trkLabArr->size(); + for (unsigned int iTrack{0}; iTrack < trkLabArr->size(); ++iTrack) { + auto lab = trkLabArr->at(iTrack); + if (!lab.isSet()) { + unaccounted++; + continue; + } + int trackID, evID, srcID; + bool fake; + lab.get(trackID, evID, srcID, fake); + if (evID < 0 || evID >= (int)info.size()) { + unaccounted++; + continue; + } + if (trackID < 0 || trackID >= (int)info[evID].size()) { + unaccounted++; + continue; + } + info[evID][trackID].isReco += !fake; + info[evID][trackID].isFake += fake; + /// We keep the best track we would keep in the data + if (recArr->at(iTrack).isBetter(info[evID][trackID].track, 1.e9)) { + info[evID][trackID].storedStatus = fake; + info[evID][trackID].track = recArr->at(iTrack); + } + + fakes += fake; + good += !fake; + } + } + std::cout << "done." << std::endl; + + std::cout << "** Some statistics:" << std::endl; + std::cout << "\t- Total number of tracks: " << total << std::endl; + std::cout << "\t- Total number of tracks not corresponding to particles: " << unaccounted << " (" << unaccounted * 100. / total << "%)" << std::endl; + std::cout << "\t- Total number of fakes: " << fakes << " (" << fakes * 100. / total << "%)" << std::endl; + std::cout << "\t- Total number of good: " << good << " (" << good * 100. / total << "%)" << std::endl; + + int nb = 100; + double xbins[nb + 1], ptcutl = 0.01, ptcuth = 10.; + double a = std::log(ptcuth / ptcutl) / nb; + for (int i = 0; i <= nb; i++) + xbins[i] = ptcutl * std::exp(i * a); + TH1D* num = new TH1D("num", ";#it{p}_{T} (GeV/#it{c});Efficiency (fake-track rate)", nb, xbins); + num->Sumw2(); + TH1D* numChi2 = new TH1D("numChi2", ";#it{p}_{T} (GeV/#it{c});Efficiency (fake-track rate)", 200, 0, 100); + + TH1D* fak = new TH1D("fak", ";#it{p}_{T} (GeV/#it{c});Fak", nb, xbins); + fak->Sumw2(); + TH1D* multiFak = new TH1D("multiFak", ";#it{p}_{T} (GeV/#it{c});Fak", nb, xbins); + multiFak->Sumw2(); + TH1D* fakChi2 = new TH1D("fakChi2", ";#it{p}_{T} (GeV/#it{c});Fak", 200, 0, 100); + + TH1D* clone = new TH1D("clone", ";#it{p}_{T} (GeV/#it{c});Clone", nb, xbins); + clone->Sumw2(); + + TH1D* den = new TH1D("den", ";#it{p}_{T} (GeV/#it{c});Den", nb, xbins); + den->Sumw2(); + + for (auto& evInfo : info) { + for (auto& part : evInfo) { + if (part.clusters != 0x7f) { + // part.clusters != 0x3f && part.clusters != 0x3f << 1 && + // part.clusters != 0x1f && part.clusters != 0x1f << 1 && part.clusters != 0x1f << 2 && + // part.clusters != 0x0f && part.clusters != 0x0f << 1 && part.clusters != 0x0f << 2 && part.clusters != 0x0f << 3) { + continue; + } + if (std::abs(part.eta) > 1.1 && !part.isPrimary) { + continue; + } + den->Fill(part.pt); + if (part.isReco) { + num->Fill(part.pt); + if (part.isReco > 1) { + for (int _i{0}; _i < part.isReco - 1; ++_i) { + clone->Fill(part.pt); + } + } + } + if (part.isFake) { + fak->Fill(part.pt); + if (part.isFake > 1) { + for (int _i{0}; _i < part.isFake - 1; ++_i) { + multiFak->Fill(part.pt); + } + } + } + } + } + + TCanvas* c1 = new TCanvas; + c1->SetLogx(); + c1->SetGridx(); + c1->SetGridy(); + TH1* sum = (TH1*)num->Clone("sum"); + sum->Add(fak); + sum->Divide(sum, den, 1, 1); + sum->SetLineColor(kBlack); + sum->Draw("hist"); + num->Divide(num, den, 1, 1, "b"); + num->Draw("histesame"); + fak->Divide(fak, den, 1, 1, "b"); + fak->SetLineColor(2); + fak->Draw("histesame"); + multiFak->Divide(multiFak, den, 1, 1, "b"); + multiFak->SetLineColor(kRed + 1); + multiFak->Draw("histsame"); + clone->Divide(clone, den, 1, 1, "b"); + clone->SetLineColor(3); + clone->Draw("histesame"); + + std::cout << "** Streaming output TTree to file ... " << std::flush; + TFile file("CheckTracksCA.root", "recreate"); + TTree tree("ParticleInfo", "ParticleInfo"); + ParticleInfo pInfo; + tree.Branch("particle", &pInfo); + 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; + } + pInfo = part; + tree.Fill(); + } + } + tree.Write(); + sum->Write("total"); + fak->Write("singleFake"); + num->Write("efficiency"); + multiFak->Write("multiFake"); + clone->Write("clones"); + file.Close(); + std::cout << " done." << std::endl; +} diff --git a/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt b/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt index c47927c93702d..e2f38e9f17244 100644 --- a/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt @@ -14,6 +14,7 @@ o2_add_library(ITStracking SOURCES src/ClusterLines.cxx src/Cluster.cxx src/ROframe.cxx + src/TimeFrame.cxx src/IOUtils.cxx src/Label.cxx src/PrimaryVertexContext.cxx diff --git a/Detectors/ITSMFT/ITS/tracking/cuda/CMakeLists.txt b/Detectors/ITSMFT/ITS/tracking/cuda/CMakeLists.txt index c47aa1f79e32b..f26c3b78dd666 100644 --- a/Detectors/ITSMFT/ITS/tracking/cuda/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/tracking/cuda/CMakeLists.txt @@ -28,4 +28,4 @@ o2_add_library(ITStrackingCUDA set_property(TARGET ${targetName} PROPERTY CUDA_SEPARABLE_COMPILATION ON) target_compile_definitions( - ${targetName} PRIVATE $) + ${targetName} PRIVATE GPUCA_O2_LIB $) diff --git a/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/ClusterLinesGPU.h b/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/ClusterLinesGPU.h index 8fa1c5a239054..9798901597cb0 100644 --- a/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/ClusterLinesGPU.h +++ b/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/ClusterLinesGPU.h @@ -17,6 +17,7 @@ #define O2_ITSMFT_TRACKING_CLUSTERLINESGPU_H_ #include "GPUCommonDef.h" +#include /// Required to properly compile MathUtils #include "ITStracking/ClusterLines.h" namespace o2 diff --git a/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/DeviceStoreNV.h b/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/DeviceStoreNV.h index 3cbc767985850..4f728f64bf633 100644 --- a/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/DeviceStoreNV.h +++ b/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/DeviceStoreNV.h @@ -16,6 +16,11 @@ #ifndef TRACKINGITSU_INCLUDE_DEVICESTORENV_H_ #define TRACKINGITSU_INCLUDE_DEVICESTORENV_H_ +#ifndef GPUCA_GPUCODE_GENRTC +#include +#include +#endif + #include "ITStracking/Cell.h" #include "ITStracking/Configuration.h" #include "ITStracking/Cluster.h" diff --git a/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/DeviceStoreVertexerGPU.h b/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/DeviceStoreVertexerGPU.h index f6214d682317d..b8a63b65f4943 100644 --- a/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/DeviceStoreVertexerGPU.h +++ b/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/DeviceStoreVertexerGPU.h @@ -19,6 +19,7 @@ #ifndef GPUCA_GPUCODE_GENRTC #include +#include #endif #include "ITStracking/Cluster.h" diff --git a/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/PrimaryVertexContextNV.h b/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/PrimaryVertexContextNV.h index 01b8b9c8d85c5..9fd1af8f714c2 100644 --- a/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/PrimaryVertexContextNV.h +++ b/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/PrimaryVertexContextNV.h @@ -18,6 +18,11 @@ #include +#ifndef GPUCA_GPUCODE_GENRTC +#include +#include +#endif + #include "ITStracking/Configuration.h" #include "ITStracking/Constants.h" #include "ITStracking/Definitions.h" diff --git a/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/TrackerTraitsNV.h b/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/TrackerTraitsNV.h index 8f0c8306cd21c..c6885a3939032 100644 --- a/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/TrackerTraitsNV.h +++ b/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/TrackerTraitsNV.h @@ -16,6 +16,10 @@ #ifndef TRACKINGITSU_INCLUDE_TRACKERTRAITSNV_H_ #define TRACKINGITSU_INCLUDE_TRACKERTRAITSNV_H_ +#ifndef GPUCA_GPUCODE_GENRTC +#include +#include +#endif #include "ITStracking/Configuration.h" #include "ITStracking/Definitions.h" #include "ITStracking/TrackerTraits.h" @@ -30,8 +34,8 @@ class PrimaryVertexContext; class TrackerTraitsNV : public TrackerTraits { public: - TrackerTraitsNV(); - ~TrackerTraitsNV() override; + TrackerTraitsNV() = default; + ~TrackerTraitsNV() override = default; void computeLayerCells() final; void computeLayerTracklets() final; diff --git a/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/Utils.h b/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/Utils.h index 6b9028bbc5947..5f3460502cb97 100644 --- a/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/Utils.h +++ b/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/Utils.h @@ -16,6 +16,11 @@ #ifndef TRACKINGITSU_INCLUDE_GPU_UTILS_H_ #define TRACKINGITSU_INCLUDE_GPU_UTILS_H_ +#ifndef GPUCA_GPUCODE_GENRTC +#include +#include +#endif + #include "GPUCommonDef.h" #include "ITStrackingCUDA/Stream.h" diff --git a/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/VertexerTraitsGPU.h b/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/VertexerTraitsGPU.h index f77339d9170d2..f9c8140ff0ae0 100644 --- a/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/VertexerTraitsGPU.h +++ b/Detectors/ITSMFT/ITS/tracking/cuda/include/ITStrackingCUDA/VertexerTraitsGPU.h @@ -71,7 +71,7 @@ inline GPUd() const int2 VertexerTraitsGPU::getBinsPhiRectWindow(const Cluster& { // This function returns the lowest PhiBin and the number of phi bins to be spanned, In the form int2{phiBinLow, PhiBinSpan} const int phiBinMin{constants::its2::getPhiBinIndex( - math_utils::getNormalizedPhiCoordinate(currentCluster.phiCoordinate - phiCut))}; + math_utils::getNormalizedPhi(currentCluster.phi - phiCut))}; const int phiBinSpan{static_cast(MATH_CEIL(phiCut * InversePhiBinSize))}; return int2{phiBinMin, phiBinSpan}; } diff --git a/Detectors/ITSMFT/ITS/tracking/cuda/src/TrackerTraitsNV.cu b/Detectors/ITSMFT/ITS/tracking/cuda/src/TrackerTraitsNV.cu index 389e42fe05726..03dc8e7bc05e5 100644 --- a/Detectors/ITSMFT/ITS/tracking/cuda/src/TrackerTraitsNV.cu +++ b/Detectors/ITSMFT/ITS/tracking/cuda/src/TrackerTraitsNV.cu @@ -13,6 +13,8 @@ /// \brief /// +#include "GPUCommonRtypes.h" + #include "ITStrackingCUDA/TrackerTraitsNV.h" #include @@ -44,9 +46,9 @@ GPU_DEVICE const int4 getBinsRect(const Cluster& currentCluster, const int layer const float z1, const float z2, float maxdeltaz, float maxdeltaphi) { const float zRangeMin = o2::gpu::GPUCommonMath::Min(z1, z2) - maxdeltaz; - const float phiRangeMin = currentCluster.phiCoordinate - maxdeltaphi; + const float phiRangeMin = currentCluster.phi - maxdeltaphi; const float zRangeMax = o2::gpu::GPUCommonMath::Max(z1, z2) + maxdeltaz; - const float phiRangeMax = currentCluster.phiCoordinate + maxdeltaphi; + const float phiRangeMax = currentCluster.phi + maxdeltaphi; if (zRangeMax < -LayersZCoordinate()[layerIndex + 1] || zRangeMin > LayersZCoordinate()[layerIndex + 1] || zRangeMin > zRangeMax) { @@ -101,9 +103,9 @@ __device__ void computeLayerTracklets(DeviceStoreNV& devStore, const int layerIn continue; }*/ - const float tanLambda{(currentCluster.zCoordinate - devStore.getPrimaryVertex().z) / currentCluster.rCoordinate}; - const float zAtRmin{tanLambda * (devStore.getRmin(layerIndex + 1) - currentCluster.rCoordinate) + currentCluster.zCoordinate}; - const float zAtRmax{tanLambda * (devStore.getRmax(layerIndex + 1) - currentCluster.rCoordinate) + currentCluster.zCoordinate}; + const float tanLambda{(currentCluster.zCoordinate - devStore.getPrimaryVertex().z) / currentCluster.radius}; + const float zAtRmin{tanLambda * (devStore.getRmin(layerIndex + 1) - currentCluster.radius) + currentCluster.zCoordinate}; + const float zAtRmax{tanLambda * (devStore.getRmax(layerIndex + 1) - currentCluster.radius) + currentCluster.zCoordinate}; const int4 selectedBinsRect{getBinsRect(currentCluster, layerIndex, zAtRmin, zAtRmax, kTrkPar.TrackletMaxDeltaZ[layerIndex], kTrkPar.TrackletMaxDeltaPhi)}; @@ -131,8 +133,8 @@ __device__ void computeLayerTracklets(DeviceStoreNV& devStore, const int layerIn const Cluster& nextCluster{nextLayerClusters[iNextLayerCluster]}; const float deltaZ{o2::gpu::GPUCommonMath::Abs( - tanLambda * (nextCluster.rCoordinate - currentCluster.rCoordinate) + currentCluster.zCoordinate - nextCluster.zCoordinate)}; - const float deltaPhi{o2::gpu::GPUCommonMath::Abs(currentCluster.phiCoordinate - nextCluster.phiCoordinate)}; + tanLambda * (nextCluster.radius - currentCluster.radius) + currentCluster.zCoordinate - nextCluster.zCoordinate)}; + const float deltaPhi{o2::gpu::GPUCommonMath::Abs(currentCluster.phi - nextCluster.phi)}; if (deltaZ < kTrkPar.TrackletMaxDeltaZ[layerIndex] && (deltaPhi < kTrkPar.TrackletMaxDeltaPhi || o2::gpu::GPUCommonMath::Abs(deltaPhi - constants::math::TwoPi) < kTrkPar.TrackletMaxDeltaPhi)) { @@ -181,8 +183,8 @@ __device__ void computeLayerCells(DeviceStoreNV& devStore, const int layerIndex, devStore.getClusters()[layerIndex][currentTracklet.firstClusterIndex]}; const Cluster& secondCellCluster{ devStore.getClusters()[layerIndex + 1][currentTracklet.secondClusterIndex]}; - const float firstCellClusterQuadraticRCoordinate{firstCellCluster.rCoordinate * firstCellCluster.rCoordinate}; - const float secondCellClusterQuadraticRCoordinate{secondCellCluster.rCoordinate * secondCellCluster.rCoordinate}; + const float firstCellClusterQuadraticRCoordinate{firstCellCluster.radius * firstCellCluster.radius}; + const float secondCellClusterQuadraticRCoordinate{secondCellCluster.radius * secondCellCluster.radius}; const float3 firstDeltaVector{secondCellCluster.xCoordinate - firstCellCluster.xCoordinate, secondCellCluster.yCoordinate - firstCellCluster.yCoordinate, secondCellClusterQuadraticRCoordinate - firstCellClusterQuadraticRCoordinate}; @@ -191,12 +193,12 @@ __device__ void computeLayerCells(DeviceStoreNV& devStore, const int layerIndex, const Tracklet& nextTracklet{devStore.getTracklets()[layerIndex + 1][iNextLayerTracklet]}; const float deltaTanLambda{o2::gpu::GPUCommonMath::Abs(currentTracklet.tanLambda - nextTracklet.tanLambda)}; - const float deltaPhi{o2::gpu::GPUCommonMath::Abs(currentTracklet.phiCoordinate - nextTracklet.phiCoordinate)}; + const float deltaPhi{o2::gpu::GPUCommonMath::Abs(currentTracklet.phi - nextTracklet.phi)}; if (deltaTanLambda < kTrkPar.CellMaxDeltaTanLambda && (deltaPhi < kTrkPar.CellMaxDeltaPhi || o2::gpu::GPUCommonMath::Abs(deltaPhi - constants::math::TwoPi) < kTrkPar.CellMaxDeltaPhi)) { const float averageTanLambda{0.5f * (currentTracklet.tanLambda + nextTracklet.tanLambda)}; - const float directionZIntersection{-averageTanLambda * firstCellCluster.rCoordinate + firstCellCluster.zCoordinate}; + const float directionZIntersection{-averageTanLambda * firstCellCluster.radius + firstCellCluster.zCoordinate}; const float deltaZ{o2::gpu::GPUCommonMath::Abs(directionZIntersection - primaryVertex.z)}; if (deltaZ < kTrkPar.CellMaxDeltaZ[layerIndex]) { @@ -204,7 +206,7 @@ __device__ void computeLayerCells(DeviceStoreNV& devStore, const int layerIndex, const Cluster& thirdCellCluster{ devStore.getClusters()[layerIndex + 2][nextTracklet.secondClusterIndex]}; - const float thirdCellClusterQuadraticRCoordinate{thirdCellCluster.rCoordinate * thirdCellCluster.rCoordinate}; + const float thirdCellClusterQuadraticRCoordinate{thirdCellCluster.radius * thirdCellCluster.radius}; const float3 secondDeltaVector{thirdCellCluster.xCoordinate - firstCellCluster.xCoordinate, thirdCellCluster.yCoordinate - firstCellCluster.yCoordinate, thirdCellClusterQuadraticRCoordinate - firstCellClusterQuadraticRCoordinate}; @@ -311,19 +313,9 @@ TrackerTraits* createTrackerTraitsNV() return new TrackerTraitsNV; } -TrackerTraitsNV::TrackerTraitsNV() -{ - mPrimaryVertexContext = new PrimaryVertexContextNV; -} - -TrackerTraitsNV::~TrackerTraitsNV() -{ - delete mPrimaryVertexContext; -} - void TrackerTraitsNV::computeLayerTracklets() { - PrimaryVertexContextNV* primaryVertexContext = static_cast(mPrimaryVertexContext); + PrimaryVertexContextNV* primaryVertexContext = static_cast(nullptr); //TODO: FIX THIS with Time Frames // cudaMemcpyToSymbol(gpu::kTrkPar, &mTrkParams, sizeof(TrackingParameters)); std::array tempSize; @@ -413,7 +405,7 @@ void TrackerTraitsNV::computeLayerTracklets() void TrackerTraitsNV::computeLayerCells() { - PrimaryVertexContextNV* primaryVertexContext = static_cast(mPrimaryVertexContext); + PrimaryVertexContextNV* primaryVertexContext = static_cast(nullptr); //TODO: FIX THIS with Time Frames std::array tempSize; std::array trackletsNum; std::array cellsNum; @@ -528,7 +520,7 @@ void TrackerTraitsNV::computeLayerCells() void TrackerTraitsNV::refitTracks(const std::vector>& tf, std::vector& tracks) { - PrimaryVertexContextNV* pvctx = static_cast(mPrimaryVertexContext); + PrimaryVertexContextNV* pvctx = static_cast(nullptr); //TODO: FIX THIS with Time Frames std::array cells; for (int iLayer = 0; iLayer < 5; iLayer++) { diff --git a/Detectors/ITSMFT/ITS/tracking/cuda/src/VertexerTraitsGPU.cu b/Detectors/ITSMFT/ITS/tracking/cuda/src/VertexerTraitsGPU.cu index c1bf8f49f5fb0..2ccc50c6a48df 100644 --- a/Detectors/ITSMFT/ITS/tracking/cuda/src/VertexerTraitsGPU.cu +++ b/Detectors/ITSMFT/ITS/tracking/cuda/src/VertexerTraitsGPU.cu @@ -40,16 +40,16 @@ namespace its using constants::its::VertexerHistogramVolume; using constants::math::TwoPi; -using math_utils::getNormalizedPhiCoordinate; +using math_utils::getNormalizedPhi; using namespace constants::its2; GPU_DEVICE const int4 getBinsRect(const Cluster& currentCluster, const int layerIndex, const float z1, float maxdeltaz, float maxdeltaphi) { const float zRangeMin = z1 - maxdeltaz; - const float phiRangeMin = currentCluster.phiCoordinate - maxdeltaphi; + const float phiRangeMin = currentCluster.phi - maxdeltaphi; const float zRangeMax = z1 + maxdeltaz; - const float phiRangeMax = currentCluster.phiCoordinate + maxdeltaphi; + const float phiRangeMax = currentCluster.phi + maxdeltaphi; if (zRangeMax < -LayersZCoordinate()[layerIndex + 1] || zRangeMin > LayersZCoordinate()[layerIndex + 1] || zRangeMin > zRangeMax) { @@ -179,7 +179,7 @@ GPUg() void trackleterKernel( const int maxRowClusterIndex{store.getIndexTable(adjacentLayerIndex)[firstBinIndex + selectedBinsRect.z - selectedBinsRect.x + 1]}; for (size_t iAdjacentCluster{(size_t)firstRowClusterIndex}; iAdjacentCluster < (size_t)maxRowClusterIndex && iAdjacentCluster < nClustersAdjacentLayer; ++iAdjacentCluster) { const Cluster& adjacentCluster = store.getClusters()[static_cast(adjacentLayerIndex)][iAdjacentCluster]; // assign-constructor may be a problem, check - if (o2::gpu::GPUCommonMath::Abs(currentCluster.phiCoordinate - adjacentCluster.phiCoordinate) < phiCut) { + if (o2::gpu::GPUCommonMath::Abs(currentCluster.phi - adjacentCluster.phi) < phiCut) { if (storedTracklets < store.getConfig().maxTrackletsPerCluster) { if (layerOrder == TrackletingLayerOrder::fromInnermostToMiddleLayer) { store.getDuplets01().emplace(stride + storedTracklets, iAdjacentCluster, currentClusterIndex, adjacentCluster, currentCluster); @@ -212,7 +212,7 @@ GPUg() void trackletSelectionKernel( for (int iTracklet12{0}; iTracklet12 < store.getNFoundTracklets(TrackletingLayerOrder::fromMiddleToOuterLayer)[currentClusterIndex]; ++iTracklet12) { for (int iTracklet01{0}; iTracklet01 < store.getNFoundTracklets(TrackletingLayerOrder::fromInnermostToMiddleLayer)[currentClusterIndex] && validTracklets < store.getConfig().maxTrackletsPerCluster; ++iTracklet01) { const float deltaTanLambda{o2::gpu::GPUCommonMath::Abs(store.getDuplets01()[stride + iTracklet01].tanLambda - store.getDuplets12()[stride + iTracklet12].tanLambda)}; - const float deltaPhi{o2::gpu::GPUCommonMath::Abs(store.getDuplets01()[stride + iTracklet01].phiCoordinate - store.getDuplets12()[stride + iTracklet12].phiCoordinate)}; + const float deltaPhi{o2::gpu::GPUCommonMath::Abs(store.getDuplets01()[stride + iTracklet01].phi - store.getDuplets12()[stride + iTracklet12].phi)}; if (deltaTanLambda < tanLambdaCut && deltaPhi < phiCut && validTracklets != store.getConfig().maxTrackletsPerCluster) { assert(store.getDuplets01()[stride + iTracklet01].secondClusterIndex == store.getDuplets12()[stride + iTracklet12].firstClusterIndex); if (!isInitRun) { diff --git a/Detectors/ITSMFT/ITS/tracking/hip/include/ITStrackingHIP/VertexerTraitsHIP.h b/Detectors/ITSMFT/ITS/tracking/hip/include/ITStrackingHIP/VertexerTraitsHIP.h index 71cef1be4c075..bd921c69ada6a 100644 --- a/Detectors/ITSMFT/ITS/tracking/hip/include/ITStrackingHIP/VertexerTraitsHIP.h +++ b/Detectors/ITSMFT/ITS/tracking/hip/include/ITStrackingHIP/VertexerTraitsHIP.h @@ -74,7 +74,7 @@ GPUdi() const int2 VertexerTraitsHIP::getBinsPhiRectWindow(const Cluster& curren { // This function returns the lowest PhiBin and the number of phi bins to be spanned, In the form int2{phiBinLow, PhiBinSpan} const int phiBinMin{constants::its2::getPhiBinIndex( - math_utils::getNormalizedPhiCoordinate(currentCluster.phiCoordinate - phiCut))}; + math_utils::getNormalizedPhi(currentCluster.phi - phiCut))}; const int phiBinSpan{static_cast(MATH_CEIL(phiCut * InversePhiBinSize))}; return int2{phiBinMin, phiBinSpan}; } diff --git a/Detectors/ITSMFT/ITS/tracking/hip/src/VertexerTraitsHIP.hip.cxx b/Detectors/ITSMFT/ITS/tracking/hip/src/VertexerTraitsHIP.hip.cxx index de83c0432c277..db95b63015376 100644 --- a/Detectors/ITSMFT/ITS/tracking/hip/src/VertexerTraitsHIP.hip.cxx +++ b/Detectors/ITSMFT/ITS/tracking/hip/src/VertexerTraitsHIP.hip.cxx @@ -43,15 +43,15 @@ using constants::its2::LayersZCoordinate; using constants::its2::PhiBins; using constants::its2::ZBins; using constants::math::TwoPi; -using math_utils::getNormalizedPhiCoordinate; +using math_utils::getNormalizedPhi; GPUhd() const int4 getBinsRect(const Cluster& currentCluster, const int layerIndex, const float z1, float maxdeltaz, float maxdeltaphi) { const float zRangeMin = z1 - maxdeltaz; - const float phiRangeMin = currentCluster.phiCoordinate - maxdeltaphi; + const float phiRangeMin = currentCluster.phi - maxdeltaphi; const float zRangeMax = z1 + maxdeltaz; - const float phiRangeMax = currentCluster.phiCoordinate + maxdeltaphi; + const float phiRangeMax = currentCluster.phi + maxdeltaphi; if (zRangeMax < -LayersZCoordinate()[layerIndex + 1] || zRangeMin > LayersZCoordinate()[layerIndex + 1] || zRangeMin > zRangeMax) { @@ -177,7 +177,7 @@ GPUg() void trackleterKernel( const int maxRowClusterIndex{store->getIndexTable(adjacentLayerIndex)[firstBinIndex + selectedBinsRect.z - selectedBinsRect.x + 1]}; for (int iAdjacentCluster{firstRowClusterIndex}; iAdjacentCluster < maxRowClusterIndex && iAdjacentCluster < nClustersAdjacentLayer; ++iAdjacentCluster) { const Cluster& adjacentCluster = store->getClusters()[static_cast(adjacentLayerIndex)][iAdjacentCluster]; // assign-constructor may be a problem, check - if (o2::gpu::GPUCommonMath::Abs(currentCluster.phiCoordinate - adjacentCluster.phiCoordinate) < phiCut) { + if (o2::gpu::GPUCommonMath::Abs(currentCluster.phi - adjacentCluster.phi) < phiCut) { if (storedTracklets < store->getConfig().maxTrackletsPerCluster) { if (layerOrder == TrackletingLayerOrder::fromInnermostToMiddleLayer) { store->getDuplets01().emplace(stride + storedTracklets, iAdjacentCluster, currentClusterIndex, adjacentCluster, currentCluster); @@ -210,7 +210,7 @@ GPUg() void trackletSelectionKernel( for (int iTracklet12{0}; iTracklet12 < store->getNFoundTracklets(TrackletingLayerOrder::fromMiddleToOuterLayer)[currentClusterIndex]; ++iTracklet12) { for (int iTracklet01{0}; iTracklet01 < store->getNFoundTracklets(TrackletingLayerOrder::fromInnermostToMiddleLayer)[currentClusterIndex] && validTracklets < store->getConfig().maxTrackletsPerCluster; ++iTracklet01) { const float deltaTanLambda{o2::gpu::GPUCommonMath::Abs(store->getDuplets01()[stride + iTracklet01].tanLambda - store->getDuplets12()[stride + iTracklet12].tanLambda)}; - const float deltaPhi{o2::gpu::GPUCommonMath::Abs(store->getDuplets01()[stride + iTracklet01].phiCoordinate - store->getDuplets12()[stride + iTracklet12].phiCoordinate)}; + const float deltaPhi{o2::gpu::GPUCommonMath::Abs(store->getDuplets01()[stride + iTracklet01].phi - store->getDuplets12()[stride + iTracklet12].phi)}; if (deltaTanLambda < tanLambdaCut && deltaPhi < phiCut && validTracklets != store->getConfig().maxTrackletsPerCluster) { assert(store->getDuplets01()[stride + iTracklet01].secondClusterIndex == store->getDuplets12()[stride + iTracklet12].firstClusterIndex); if (!isInitRun) { diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Cluster.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Cluster.h index 9b97bd2d7cae8..81b8ccb37dcfe 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Cluster.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Cluster.h @@ -41,8 +41,8 @@ struct Cluster final { float xCoordinate; // = -999.f; float yCoordinate; // = -999.f; float zCoordinate; // = -999.f; - float phiCoordinate; // = -999.f; - float rCoordinate; // = -999.f; + float phi; // = -999.f; + float radius; // = -999.f; int clusterId; // = -1; int indexTableBinIndex; // = -1; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h index d3fba50eb4b8c..5ff3cb0225785 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h @@ -48,12 +48,28 @@ class Configuration : public Param struct TrackingParameters { TrackingParameters& operator=(const TrackingParameters& t) = default; + void CopyCuts(TrackingParameters& other, float scale = 1.) + { + TrackletMaxDeltaPhi = other.TrackletMaxDeltaPhi * scale; + for (unsigned int ii{0}; ii < TrackletMaxDeltaZ.size(); ++ii) { + TrackletMaxDeltaZ[ii] = other.TrackletMaxDeltaZ[ii] * scale; + } + CellMaxDeltaTanLambda = other.CellMaxDeltaTanLambda * scale; + for (unsigned int ii{0}; ii < CellMaxDCA.size(); ++ii) { + CellMaxDCA[ii] = other.CellMaxDCA[ii] * scale; + } + for (unsigned int ii{0}; ii < NeighbourMaxDeltaCurvature.size(); ++ii) { + NeighbourMaxDeltaCurvature[ii] = other.NeighbourMaxDeltaCurvature[ii] * scale; + NeighbourMaxDeltaN[ii] = other.NeighbourMaxDeltaN[ii] * scale; + } + } int CellMinimumLevel(); int CellsPerRoad() const { return NLayers - 2; } int TrackletsPerRoad() const { return NLayers - 1; } int NLayers = 7; + int DeltaROF = 1; std::vector LayerZ = {16.333f + 1, 16.333f + 1, 16.333f + 1, 42.140f + 1, 42.140f + 1, 73.745f + 1, 73.745f + 1}; std::vector LayerRadii = {2.33959f, 3.14076f, 3.91924f, 19.6213f, 24.5597f, 34.388f, 39.3329f}; int ZBins{256}; @@ -75,7 +91,7 @@ struct TrackingParameters { std::vector NeighbourMaxDeltaN = {0.002f, 0.0090f, 0.002f, 0.005f}; /// Fitter parameters bool UseMatBudLUT = false; - std::array FitIterationMaxChi2 = {o2::constants::math::VeryBig, o2::constants::math::VeryBig}; + std::array FitIterationMaxChi2 = {100, 50}; }; struct MemoryParameters { diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Constants.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Constants.h index 1f045f3aa9d13..25825dace2404 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Constants.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Constants.h @@ -113,6 +113,11 @@ namespace pdgcodes constexpr int PionCode{211}; } } // namespace constants +// typedef std::array, +// constants::its::TrackletsPerRoad> index_table_t; +#ifndef __OPENCL__ /// FIXME: this is for compatibility with OCL +typedef std::vector> index_table_t; +#endif } // namespace its } // namespace o2 diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Label.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Label.h index bb39c9349e273..ec45e6587a974 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Label.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Label.h @@ -28,7 +28,7 @@ struct Label final { int monteCarloId; float transverseMomentum; - float phiCoordinate; + float phi; float pseudorapidity; int pdgCode; int numberOfClusters; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h index 16b839afba49e..b9a2b599f0ba8 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h @@ -35,9 +35,9 @@ namespace its namespace math_utils { -GPUhdni() float calculatePhiCoordinate(const float, const float); -GPUhdni() float calculateRCoordinate(const float, const float); -GPUhdni() constexpr float getNormalizedPhiCoordinate(const float); +GPUhdni() float computePhi(const float, const float); +GPUhdni() float hypot(const float, const float); +GPUhdni() constexpr float getNormalizedPhi(const float); GPUhdni() constexpr float3 crossProduct(const float3&, const float3&); GPUhdni() float computeCurvature(float x1, float y1, float x2, float y2, float x3, float y3); GPUhdni() float computeCurvatureCentreX(float x1, float y1, float x2, float y2, float x3, float y3); @@ -45,22 +45,21 @@ GPUhdni() float computeTanDipAngle(float x1, float y1, float x2, float y2, float } // namespace math_utils -GPUhdi() float math_utils::calculatePhiCoordinate(const float xCoordinate, const float yCoordinate) +GPUhdi() float math_utils::computePhi(const float x, const float y) { //return o2::gpu::CAMath::ATan2(-yCoordinate, -xCoordinate) + constants::math::Pi; - return o2::math_utils::fastATan2(-yCoordinate, -xCoordinate) + constants::math::Pi; + return o2::math_utils::fastATan2(-y, -x) + constants::math::Pi; } -GPUhdi() float math_utils::calculateRCoordinate(const float xCoordinate, const float yCoordinate) +GPUhdi() float math_utils::hypot(const float x, const float y) { - return o2::gpu::CAMath::Sqrt(xCoordinate * xCoordinate + yCoordinate * yCoordinate); + return o2::gpu::CAMath::Sqrt(x * x + y * y); } -GPUhdi() constexpr float math_utils::getNormalizedPhiCoordinate(const float phiCoordinate) +GPUhdi() constexpr float math_utils::getNormalizedPhi(const float phi) { - return (phiCoordinate < 0) - ? phiCoordinate + constants::math::TwoPi - : (phiCoordinate > constants::math::TwoPi) ? phiCoordinate - constants::math::TwoPi : phiCoordinate; + return (phi < 0) ? phi + constants::math::TwoPi : (phi > constants::math::TwoPi) ? phi - constants::math::TwoPi + : phi; } GPUhdi() constexpr float3 math_utils::crossProduct(const float3& firstVector, const float3& secondVector) diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h new file mode 100644 index 0000000000000..40f723b895f0e --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h @@ -0,0 +1,368 @@ +// 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. +// +// 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 TimeFrame.h +/// \brief +/// + +#ifndef TRACKINGITSU_INCLUDE_TIMEFRAME_H_ +#define TRACKINGITSU_INCLUDE_TIMEFRAME_H_ + +#include +#include +#include +#include +#include + +#include "DataFormatsITS/TrackITS.h" + +#include "ITStracking/Cell.h" +#include "ITStracking/Cluster.h" +#include "ITStracking/Configuration.h" +#include "ITStracking/Constants.h" +#include "ITStracking/Definitions.h" +#include "ITStracking/Road.h" +#include "ITStracking/Tracklet.h" +#include "ITStracking/IndexTableUtils.h" + +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" + +#include "ReconstructionDataFormats/Vertex.h" + +namespace o2 +{ + +namespace itsmft +{ +class Cluster; +class CompClusterExt; +class TopologyDictionary; +class ROFRecord; +} // namespace itsmft + +namespace its +{ + +using Vertex = o2::dataformats::Vertex>; + +class TimeFrame final +{ + public: + TimeFrame(int nLayers = 7); + const Vertex& getPrimaryVertex(const int) const; + gsl::span getPrimaryVertices(int tf) const; + gsl::span getPrimaryVertices(int romin, int romax) const; + int getPrimaryVerticesNum(int rofID = -1) const; + void addPrimaryVertices(const std::vector& vertices); + int loadROFrameData(const o2::itsmft::ROFRecord& rof, gsl::span clusters, + const dataformats::MCTruthContainer* mcLabels = nullptr); + + int loadROFrameData(gsl::span rofs, gsl::span clusters, gsl::span::iterator& pattIt, + const itsmft::TopologyDictionary& dict, const dataformats::MCTruthContainer* mcLabels = nullptr); + int getTotalClusters() const; + bool empty() const; + + int getSortedIndex(int rof, int layer, int i) const; + int getNrof() const; + float getBeamX() const; + float getBeamY() const; + + float getMinR(int layer) const { return mMinR[layer]; } + float getMaxR(int layer) const { return mMaxR[layer]; } + + gsl::span getClustersOnLayer(int rofId, int layerId); + gsl::span getClustersOnLayer(int rofId, int layerId) const; + gsl::span getUnsortedClustersOnLayer(int rofId, int layerId) const; + index_table_t& getIndexTables(int tf); + const std::vector& getTrackingFrameInfoOnLayer(int layerId) const; + + const TrackingFrameInfo& getClusterTrackingFrameInfo(int layerId, const Cluster& cl) const; + const MCCompLabel& getClusterLabels(int layerId, const Cluster& cl) const; + const MCCompLabel& getClusterLabels(int layerId, const int clId) const; + int getClusterExternalIndex(int layerId, const int clId) const; + + std::vector& getTrackletsLabel(int layer) { return mTrackletLabels[layer]; } + std::vector& getCellsLabel(int layer) { return mCellLabels[layer]; } + + bool hasMCinformation() const; + void initialise(const int iteration, const MemoryParameters& memParam, const TrackingParameters& trkParam); + + bool isClusterUsed(int layer, int clusterId) const; + void markUsedCluster(int layer, int clusterId); + + std::vector>& getTracklets(); + std::vector>& getTrackletsLookupTable(); + + std::vector>& getClusters(); + std::vector>& getUnsortedClusters(); + int getClusterROF(int iLayer, int iCluster); + std::vector>& getCells(); + std::vector>& getCellsLookupTable(); + std::vector>>& getCellsNeighbours(); + std::vector& getRoads(); + std::vector& getTracks(int rof) { return mTracks[rof]; } + std::vector& getTracksLabel(int rof) { return mTracksLabel[rof]; } + + void initialiseRoadLabels(); + void setRoadLabel(int i, const unsigned long long& lab, bool fake); + const unsigned long long& getRoadLabel(int i) const; + bool isRoadFake(int i) const; + + void clear(); + + /// Debug and printing + void checkTrackletLUTs(); + void printROFoffsets(); + void printVertices(); + void printTrackletLUTonLayer(int i); + void printCellLUTonLayer(int i); + void printTrackletLUTs(); + void printCellLUTs(); + + IndexTableUtils mIndexTableUtils; + + private: + template + void addClusterToLayer(int layer, T&&... args); + template + void addTrackingFrameInfoToLayer(int layer, T&&... args); + void addClusterLabelToLayer(int layer, const MCCompLabel label); + void addClusterExternalIndexToLayer(int layer, const int idx); + + int mNrof = 0; + int mBeamPosWeight = 0; + float mBeamPos[2] = {0.f, 0.f}; + std::vector mMinR; + std::vector mMaxR; + std::vector mROframesPV = {0}; + std::vector> mROframesClusters; + std::vector mPrimaryVertices; + std::vector> mClusters; + std::vector> mUnsortedClusters; + std::vector> mUsedClusters; + std::vector> mTrackingFrameInfo; + std::vector> mClusterLabels; + std::vector> mTrackletLabels; + std::vector> mCellLabels; + std::vector> mClusterExternalIndices; + std::vector> mCells; + std::vector> mCellsLookupTable; + std::vector>> mCellsNeighbours; + std::vector mRoads; + std::vector> mTracksLabel; + std::vector> mTracks; + + std::vector mIndexTables; + std::vector> mTracklets; + std::vector> mTrackletsLookupTable; + + std::vector> mRoadLabels; +}; + +inline const Vertex& TimeFrame::getPrimaryVertex(const int vertexIndex) const { return mPrimaryVertices[vertexIndex]; } + +inline gsl::span TimeFrame::getPrimaryVertices(int tf) const +{ + const int start = mROframesPV[tf]; + const int stop = tf >= mNrof - 1 ? mNrof : tf + 1; + return {&mPrimaryVertices[start], static_cast::size_type>(mROframesPV[stop] - mROframesPV[tf])}; +} + +inline gsl::span TimeFrame::getPrimaryVertices(int romin, int romax) const +{ + return {&mPrimaryVertices[mROframesPV[romin]], static_cast::size_type>(mROframesPV[romax + 1] - mROframesPV[romin])}; +} + +inline int TimeFrame::getPrimaryVerticesNum(int rofID) const +{ + return rofID < 0 ? mPrimaryVertices.size() : mROframesPV[rofID + 1] - mROframesPV[rofID]; +} + +inline bool TimeFrame::empty() const { return getTotalClusters() == 0; } + +inline int TimeFrame::getSortedIndex(int rof, int layer, int index) const { return mROframesClusters[layer][rof] + index; } + +inline int TimeFrame::getNrof() const { return mNrof; }; + +inline float TimeFrame::getBeamX() const { return mBeamPos[0]; } + +inline float TimeFrame::getBeamY() const { return mBeamPos[1]; } + +inline gsl::span TimeFrame::getClustersOnLayer(int rofId, int layerId) +{ + if (rofId < 0 || rofId >= mNrof) { + return gsl::span(); + } + int startIdx{mROframesClusters[layerId][rofId]}; + return {&mClusters[layerId][startIdx], static_cast::size_type>(mROframesClusters[layerId][rofId + 1] - startIdx)}; +} + +inline gsl::span TimeFrame::getClustersOnLayer(int rofId, int layerId) const +{ + if (rofId < 0 || rofId >= mNrof) { + return gsl::span(); + } + int startIdx{mROframesClusters[layerId][rofId]}; + return {&mClusters[layerId][startIdx], static_cast::size_type>(mROframesClusters[layerId][rofId + 1] - startIdx)}; +} + +inline int TimeFrame::getClusterROF(int iLayer, int iCluster) +{ + return std::lower_bound(mROframesClusters[iLayer].begin(), mROframesClusters[iLayer].end(), iCluster + 1) - mROframesClusters[iLayer].begin() - 1; +} + +inline gsl::span TimeFrame::getUnsortedClustersOnLayer(int rofId, int layerId) const +{ + if (rofId < 0 || rofId >= mNrof) { + return gsl::span(); + } + int startIdx{mROframesClusters[layerId][rofId]}; + return {&mUnsortedClusters[layerId][startIdx], static_cast::size_type>(mROframesClusters[layerId][rofId + 1] - startIdx)}; +} + +inline const std::vector& TimeFrame::getTrackingFrameInfoOnLayer(int layerId) const +{ + return mTrackingFrameInfo[layerId]; +} + +inline const TrackingFrameInfo& TimeFrame::getClusterTrackingFrameInfo(int layerId, const Cluster& cl) const +{ + return mTrackingFrameInfo[layerId][cl.clusterId]; +} + +inline const MCCompLabel& TimeFrame::getClusterLabels(int layerId, const Cluster& cl) const +{ + return mClusterLabels[layerId][cl.clusterId]; +} + +inline const MCCompLabel& TimeFrame::getClusterLabels(int layerId, const int clId) const +{ + return mClusterLabels[layerId][clId]; +} + +inline int TimeFrame::getClusterExternalIndex(int layerId, const int clId) const +{ + return mClusterExternalIndices[layerId][clId]; +} + +inline index_table_t& TimeFrame::getIndexTables(int tf) +{ + return mIndexTables[tf]; +} + +template +void TimeFrame::addClusterToLayer(int layer, T&&... values) +{ + mUnsortedClusters[layer].emplace_back(std::forward(values)...); +} + +template +void TimeFrame::addTrackingFrameInfoToLayer(int layer, T&&... values) +{ + mTrackingFrameInfo[layer].emplace_back(std::forward(values)...); +} + +inline void TimeFrame::addClusterLabelToLayer(int layer, const MCCompLabel label) { mClusterLabels[layer].emplace_back(label); } + +inline void TimeFrame::addClusterExternalIndexToLayer(int layer, const int idx) +{ + mClusterExternalIndices[layer].push_back(idx); +} + +inline void TimeFrame::clear() +{ + for (unsigned int iL = 0; iL < mClusters.size(); ++iL) { + mClusters[iL].clear(); + mTrackingFrameInfo[iL].clear(); + mClusterLabels[iL].clear(); + mClusterExternalIndices[iL].clear(); + } + mPrimaryVertices.clear(); +} + +inline bool TimeFrame::hasMCinformation() const +{ + for (const auto& vect : mClusterLabels) { + if (!vect.empty()) { + return true; + } + } + return false; +} + +inline bool TimeFrame::isClusterUsed(int layer, int clusterId) const +{ + return mUsedClusters[layer][clusterId]; +} + +inline void TimeFrame::markUsedCluster(int layer, int clusterId) { mUsedClusters[layer][clusterId] = true; } + +inline std::vector>& TimeFrame::getTracklets() +{ + return mTracklets; +} + +inline std::vector>& TimeFrame::getTrackletsLookupTable() +{ + return mTrackletsLookupTable; +} + +inline void TimeFrame::initialiseRoadLabels() +{ + mRoadLabels.clear(); + mRoadLabels.resize(mRoads.size()); +} + +inline void TimeFrame::setRoadLabel(int i, const unsigned long long& lab, bool fake) +{ + mRoadLabels[i].first = lab; + mRoadLabels[i].second = fake; +} + +inline const unsigned long long& TimeFrame::getRoadLabel(int i) const +{ + return mRoadLabels[i].first; +} + +inline bool TimeFrame::isRoadFake(int i) const +{ + return mRoadLabels[i].second; +} + +inline std::vector>& TimeFrame::getClusters() +{ + return mClusters; +} + +inline std::vector>& TimeFrame::getUnsortedClusters() +{ + return mUnsortedClusters; +} + +inline std::vector>& TimeFrame::getCells() { return mCells; } + +inline std::vector>& TimeFrame::getCellsLookupTable() +{ + return mCellsLookupTable; +} + +inline std::vector>>& TimeFrame::getCellsNeighbours() +{ + return mCellsNeighbours; +} + +inline std::vector& TimeFrame::getRoads() { return mRoads; } + +} // namespace its +} // namespace o2 + +#endif /* TRACKINGITSU_INCLUDE_TimeFrame_H_ */ diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h index 56b4d5a14212c..47fdadb44a0cb 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h @@ -24,6 +24,7 @@ #include #include #include +#include #include "ITStracking/Configuration.h" #include "DetectorsBase/MatLayerCylSet.h" @@ -31,15 +32,13 @@ #include "ITStracking/Definitions.h" #include "ITStracking/ROframe.h" #include "ITStracking/MathUtils.h" -#include "ITStracking/PrimaryVertexContext.h" #include "DetectorsBase/Propagator.h" +#include "ITStracking/TimeFrame.h" #include "ITStracking/Road.h" #include "DataFormatsITS/TrackITS.h" #include "SimulationDataFormat/MCCompLabel.h" -#include "Framework/Logger.h" - #ifdef CA_DEBUG #include "ITStracking/StandaloneDebugger.h" #endif @@ -54,7 +53,7 @@ class GPUChainITS; namespace its { -class PrimaryVertexContext; +class TimeFrame; class TrackerTraits; class Tracker @@ -70,10 +69,7 @@ class Tracker void setBz(float bz); float getBz() const; - std::vector& getTracks(); - auto& getTrackLabels() { return mTrackLabels; } - - void clustersToTracks(const ROframe&, std::ostream& = std::cout); + void clustersToTracks(std::function = [](std::string s) { std::cout << s << std::endl; }); void setROFrame(std::uint32_t f) { mROFrame = f; } std::uint32_t getROFrame() const { return mROFrame; } @@ -86,23 +82,23 @@ class Tracker track::TrackParCov buildTrackSeed(const Cluster& cluster1, const Cluster& cluster2, const Cluster& cluster3, const TrackingFrameInfo& tf3); template - void initialisePrimaryVertexContext(T&&... args); + void initialiseTimeFrame(T&&... args); void computeTracklets(); void computeCells(); void findCellsNeighbours(int& iteration); void findRoads(int& iteration); - void findTracks(const ROframe& ev); - bool fitTrack(const ROframe& event, TrackITSExt& track, int start, int end, int step, const float chi2cut = o2::constants::math::VeryBig); + void findTracks(); + bool fitTrack(TrackITSExt& track, int start, int end, int step, const float chi2cut = o2::constants::math::VeryBig, const float maxQoverPt = o2::constants::math::VeryBig); void traverseCellsTree(const int, const int); - void computeRoadsMClabels(const ROframe&); - void computeTracksMClabels(const ROframe&); - void rectifyClusterIndices(const ROframe& event); + void computeRoadsMClabels(); + void computeTracksMClabels(); + void rectifyClusterIndices(); template - float evaluateTask(void (Tracker::*)(T...), const char*, std::ostream& ostream, T&&... args); + float evaluateTask(void (Tracker::*)(T...), const char*, std::function logger, T&&... args); TrackerTraits* mTraits = nullptr; /// Observer pointer, not owned by this class - PrimaryVertexContext* mPrimaryVertexContext = nullptr; /// Observer pointer, not owned by this class + TimeFrame* mTimeFrame = nullptr; /// Observer pointer, not owned by this class std::vector mMemParams; std::vector mTrkParams; @@ -111,8 +107,6 @@ class Tracker o2::base::PropagatorImpl::MatCorrType mCorrType = o2::base::PropagatorImpl::MatCorrType::USEMatCorrLUT; float mBz = 5.f; std::uint32_t mROFrame = 0; - std::vector mTracks; - std::vector mTrackLabels; o2::gpu::GPUChainITS* mRecoChain = nullptr; #ifdef CA_DEBUG @@ -137,18 +131,13 @@ inline void Tracker::setBz(float bz) } template -void Tracker::initialisePrimaryVertexContext(T&&... args) -{ - mPrimaryVertexContext->initialise(std::forward(args)...); -} - -inline std::vector& Tracker::getTracks() +void Tracker::initialiseTimeFrame(T&&... args) { - return mTracks; + mTimeFrame->initialise(std::forward(args)...); } template -float Tracker::evaluateTask(void (Tracker::*task)(T...), const char* taskName, std::ostream& ostream, +float Tracker::evaluateTask(void (Tracker::*task)(T...), const char* taskName, std::function logger, T&&... args) { float diff{0.f}; @@ -161,13 +150,13 @@ float Tracker::evaluateTask(void (Tracker::*task)(T...), const char* taskName, s std::chrono::duration diff_t{end - start}; diff = diff_t.count(); - if (fair::Logger::Logging(fair::Severity::info)) { - if (taskName == nullptr) { - ostream << diff << "\t"; - } else { - ostream << std::setw(2) << " - " << taskName << " completed in: " << diff << " ms" << std::endl; - } + std::stringstream sstream; + if (taskName == nullptr) { + sstream << diff << "\t"; + } else { + sstream << std::setw(2) << " - " << taskName << " completed in: " << diff << " ms"; } + logger(sstream.str()); } else { (this->*task)(std::forward(args)...); } diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h index 9bb03bb328a42..6fce049e70665 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h @@ -29,7 +29,7 @@ #include "ITStracking/Configuration.h" #include "ITStracking/Definitions.h" #include "ITStracking/MathUtils.h" -#include "ITStracking/PrimaryVertexContext.h" +#include "ITStracking/TimeFrame.h" #include "ITStracking/Road.h" namespace o2 @@ -63,10 +63,10 @@ class TrackerTraits virtual void refitTracks(const std::vector>&, std::vector&){}; void UpdateTrackingParameters(const TrackingParameters& trkPar); - PrimaryVertexContext* getPrimaryVertexContext() { return mPrimaryVertexContext; } + TimeFrame* getTimeFrame() { return mTimeFrame; } protected: - PrimaryVertexContext* mPrimaryVertexContext; + TimeFrame* mTimeFrame; TrackingParameters mTrkParams; o2::gpu::GPUChainITS* mChain = nullptr; @@ -82,9 +82,9 @@ inline GPU_DEVICE const int4 TrackerTraits::getBinsRect(const Cluster& currentCl const float z1, const float z2, float maxdeltaz, float maxdeltaphi) { const float zRangeMin = o2::gpu::GPUCommonMath::Min(z1, z2) - maxdeltaz; - const float phiRangeMin = currentCluster.phiCoordinate - maxdeltaphi; + const float phiRangeMin = currentCluster.phi - maxdeltaphi; const float zRangeMax = o2::gpu::GPUCommonMath::Max(z1, z2) + maxdeltaz; - const float phiRangeMax = currentCluster.phiCoordinate + maxdeltaphi; + const float phiRangeMax = currentCluster.phi + maxdeltaphi; if (zRangeMax < -mTrkParams.LayerZ[layerIndex + 1] || zRangeMin > mTrkParams.LayerZ[layerIndex + 1] || zRangeMin > zRangeMax) { @@ -92,11 +92,11 @@ inline GPU_DEVICE const int4 TrackerTraits::getBinsRect(const Cluster& currentCl return getEmptyBinsRect(); } - const IndexTableUtils& utils{mPrimaryVertexContext->mIndexTableUtils}; + const IndexTableUtils& utils{mTimeFrame->mIndexTableUtils}; return int4{o2::gpu::GPUCommonMath::Max(0, utils.getZBinIndex(layerIndex + 1, zRangeMin)), - utils.getPhiBinIndex(math_utils::getNormalizedPhiCoordinate(phiRangeMin)), + utils.getPhiBinIndex(math_utils::getNormalizedPhi(phiRangeMin)), o2::gpu::GPUCommonMath::Min(mTrkParams.ZBins - 1, utils.getZBinIndex(layerIndex + 1, zRangeMax)), - utils.getPhiBinIndex(math_utils::getNormalizedPhiCoordinate(phiRangeMax))}; + utils.getPhiBinIndex(math_utils::getNormalizedPhi(phiRangeMax))}; } } // namespace its } // namespace o2 diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraitsCPU.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraitsCPU.h index 3fbb208e06d27..c6ef2e511e250 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraitsCPU.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraitsCPU.h @@ -29,7 +29,7 @@ #include "ITStracking/Configuration.h" #include "ITStracking/Definitions.h" #include "ITStracking/MathUtils.h" -#include "ITStracking/PrimaryVertexContext.h" +#include "ITStracking/TimeFrame.h" #include "ITStracking/Road.h" namespace o2 @@ -40,8 +40,8 @@ namespace its class TrackerTraitsCPU : public TrackerTraits { public: - TrackerTraitsCPU() { mPrimaryVertexContext = new PrimaryVertexContext; } - ~TrackerTraitsCPU() override { delete mPrimaryVertexContext; } + TrackerTraitsCPU(TimeFrame* tf = nullptr) { mTimeFrame = tf; } //TODO: the TimeFrame pointer should be given by the tracker that is the external interface + ~TrackerTraitsCPU() override { delete mTimeFrame; } void computeLayerCells() final; void computeLayerTracklets() final; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracklet.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracklet.h index e8cc02a81dd79..377ca7a6906f2 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracklet.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracklet.h @@ -28,7 +28,7 @@ namespace its struct Tracklet final { Tracklet(); - GPUdi() Tracklet(const int, const int, const Cluster&, const Cluster&); + GPUdi() Tracklet(const int, const int, const Cluster&, const Cluster&, int rof0, int rof1); #ifdef _ALLOW_DEBUG_TREES_ITS_ unsigned char isEmpty() const; void dump(); @@ -38,22 +38,24 @@ struct Tracklet final { int firstClusterIndex; int secondClusterIndex; float tanLambda; - float phiCoordinate; + float phi; + unsigned short rof[2]; }; -inline Tracklet::Tracklet() : firstClusterIndex{0}, secondClusterIndex{0}, tanLambda{0.0f}, phiCoordinate{0.0f} +inline Tracklet::Tracklet() : firstClusterIndex{0}, secondClusterIndex{0}, tanLambda{0.0f}, phi{0.0f} { // Nothing to do } GPUdi() Tracklet::Tracklet(const int firstClusterOrderingIndex, const int secondClusterOrderingIndex, - const Cluster& firstCluster, const Cluster& secondCluster) + const Cluster& firstCluster, const Cluster& secondCluster, int rof0 = -1, int rof1 = -1) : firstClusterIndex{firstClusterOrderingIndex}, secondClusterIndex{secondClusterOrderingIndex}, tanLambda{(firstCluster.zCoordinate - secondCluster.zCoordinate) / - (firstCluster.rCoordinate - secondCluster.rCoordinate)}, - phiCoordinate{o2::gpu::GPUCommonMath::ATan2(firstCluster.yCoordinate - secondCluster.yCoordinate, - firstCluster.xCoordinate - secondCluster.xCoordinate)} + (firstCluster.radius - secondCluster.radius)}, + phi{o2::gpu::GPUCommonMath::ATan2(firstCluster.yCoordinate - secondCluster.yCoordinate, + firstCluster.xCoordinate - secondCluster.xCoordinate)}, + rof{static_cast(rof0), static_cast(rof1)} { // Nothing to do } @@ -61,7 +63,7 @@ GPUdi() Tracklet::Tracklet(const int firstClusterOrderingIndex, const int second #ifdef _ALLOW_DEBUG_TREES_ITS_ inline unsigned char Tracklet::isEmpty() const { - return !firstClusterIndex && !secondClusterIndex && !tanLambda && !phiCoordinate; + return !firstClusterIndex && !secondClusterIndex && !tanLambda && !phi; } inline unsigned char Tracklet::operator<(const Tracklet& t) const @@ -80,7 +82,7 @@ inline void Tracklet::dump() std::cout << "firstClusterIndex: " << firstClusterIndex << std::endl; std::cout << "secondClusterIndex: " << secondClusterIndex << std::endl; std::cout << "tanLambda: " << tanLambda << std::endl; - std::cout << "phiCoordinate: " << phiCoordinate << std::endl; + std::cout << "phi: " << phi << std::endl; } #endif diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h index 51fda44ffd8b6..061e5e692ba23 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h @@ -60,7 +60,8 @@ class Vertexer std::vector exportVertices(); VertexerTraits* getTraits() const { return mTraits; }; - float clustersToVertices(ROframe&, const bool useMc = false, std::ostream& = std::cout); + float clustersToVertices( + ROframe&, const bool useMc = false, std::function = [](std::string s) { std::cout << s << std::endl; }); void filterMCTracklets(); void validateTracklets(); @@ -77,7 +78,7 @@ class Vertexer // Utils void dumpTraits(); template - float evaluateTask(void (Vertexer::*)(T...), const char*, std::ostream& ostream, T&&... args); + float evaluateTask(void (Vertexer::*)(T...), const char*, std::function logger, T&&... args); // debug void setDebugCombinatorics(); @@ -140,9 +141,6 @@ inline std::vector Vertexer::exportVertices() { std::vector vertices; for (auto& vertex : mTraits->getVertices()) { - if (fair::Logger::Logging(fair::Severity::info)) { - std::cout << "\t\tFound vertex with: " << std::setw(6) << vertex.mContributors << " contributors" << std::endl; - } vertices.emplace_back(o2::math_utils::Point3D(vertex.mX, vertex.mY, vertex.mZ), vertex.mRMS2, vertex.mContributors, vertex.mAvgDistance2); vertices.back().setTimeStamp(vertex.mTimeStamp); } @@ -150,7 +148,7 @@ inline std::vector Vertexer::exportVertices() } template -float Vertexer::evaluateTask(void (Vertexer::*task)(T...), const char* taskName, std::ostream& ostream, +float Vertexer::evaluateTask(void (Vertexer::*task)(T...), const char* taskName, std::function logger, T&&... args) { float diff{0.f}; @@ -163,13 +161,13 @@ float Vertexer::evaluateTask(void (Vertexer::*task)(T...), const char* taskName, std::chrono::duration diff_t{end - start}; diff = diff_t.count(); - if (fair::Logger::Logging(fair::Severity::info)) { - if (taskName == nullptr) { - ostream << diff << "\t"; - } else { - ostream << std::setw(2) << " - " << taskName << " completed in: " << diff << " ms" << std::endl; - } + std::stringstream sstream; + if (taskName == nullptr) { + sstream << diff << "\t"; + } else { + sstream << std::setw(2) << " - " << taskName << " completed in: " << diff << " ms"; } + logger(sstream.str()); } else { (this->*task)(std::forward(args)...); } diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexerTraits.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexerTraits.h index 3d7da70dfc84c..ed614254eab1e 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexerTraits.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexerTraits.h @@ -210,8 +210,8 @@ GPUhdi() const int2 VertexerTraits::getPhiBins(float phi, float dPhi) GPUhdi() const int2 VertexerTraits::getPhiBins(float phi, float dPhi, const IndexTableUtils& utils) { - return int2{utils.getPhiBinIndex(math_utils::getNormalizedPhiCoordinate(phi - dPhi)), - utils.getPhiBinIndex(math_utils::getNormalizedPhiCoordinate(phi + dPhi))}; + return int2{utils.getPhiBinIndex(math_utils::getNormalizedPhi(phi - dPhi)), + utils.getPhiBinIndex(math_utils::getNormalizedPhi(phi + dPhi))}; } GPUhdi() const int4 VertexerTraits::getBinsRect(const Cluster& currentCluster, const int layerIndex, @@ -219,9 +219,9 @@ GPUhdi() const int4 VertexerTraits::getBinsRect(const Cluster& currentCluster, c const IndexTableUtils& utils) { const float zRangeMin = directionZIntersection - 2 * maxdeltaz; - const float phiRangeMin = currentCluster.phiCoordinate - maxdeltaphi; + const float phiRangeMin = currentCluster.phi - maxdeltaphi; const float zRangeMax = directionZIntersection + 2 * maxdeltaz; - const float phiRangeMax = currentCluster.phiCoordinate + maxdeltaphi; + const float phiRangeMax = currentCluster.phi + maxdeltaphi; if (zRangeMax < -utils.getLayerZ(layerIndex + 1) || zRangeMin > utils.getLayerZ(layerIndex + 1) || zRangeMin > zRangeMax) { @@ -230,9 +230,9 @@ GPUhdi() const int4 VertexerTraits::getBinsRect(const Cluster& currentCluster, c } return int4{o2::gpu::GPUCommonMath::Max(0, utils.getZBinIndex(layerIndex + 1, zRangeMin)), - utils.getPhiBinIndex(math_utils::getNormalizedPhiCoordinate(phiRangeMin)), + utils.getPhiBinIndex(math_utils::getNormalizedPhi(phiRangeMin)), o2::gpu::GPUCommonMath::Min(utils.getNzBins() - 1, utils.getZBinIndex(layerIndex + 1, zRangeMax)), - utils.getPhiBinIndex(math_utils::getNormalizedPhiCoordinate(phiRangeMax))}; + utils.getPhiBinIndex(math_utils::getNormalizedPhi(phiRangeMax))}; } GPUhdi() const int4 VertexerTraits::getBinsRect(const Cluster& currentCluster, const int layerIndex, diff --git a/Detectors/ITSMFT/ITS/tracking/src/Cluster.cxx b/Detectors/ITSMFT/ITS/tracking/src/Cluster.cxx index 9e92af4f510d9..f4bac7dc62c7e 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Cluster.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Cluster.cxx @@ -22,16 +22,16 @@ namespace o2 namespace its { -using math_utils::calculatePhiCoordinate; -using math_utils::calculateRCoordinate; -using math_utils::getNormalizedPhiCoordinate; +using math_utils::computePhi; +using math_utils::getNormalizedPhi; +using math_utils::hypot; Cluster::Cluster(const float x, const float y, const float z, const int index) : xCoordinate{x}, yCoordinate{y}, zCoordinate{z}, - phiCoordinate{getNormalizedPhiCoordinate(calculatePhiCoordinate(x, y))}, - rCoordinate{calculateRCoordinate(x, y)}, + phi{getNormalizedPhi(computePhi(x, y))}, + radius{hypot(x, y)}, clusterId{index}, indexTableBinIndex{0} { @@ -42,11 +42,11 @@ Cluster::Cluster(const int layerIndex, const IndexTableUtils& utils, const Clust : xCoordinate{other.xCoordinate}, yCoordinate{other.yCoordinate}, zCoordinate{other.zCoordinate}, - phiCoordinate{getNormalizedPhiCoordinate(calculatePhiCoordinate(other.xCoordinate, other.yCoordinate))}, - rCoordinate{calculateRCoordinate(other.xCoordinate, other.yCoordinate)}, + phi{getNormalizedPhi(computePhi(other.xCoordinate, other.yCoordinate))}, + radius{hypot(other.xCoordinate, other.yCoordinate)}, clusterId{other.clusterId}, indexTableBinIndex{utils.getBinIndex(utils.getZBinIndex(layerIndex, zCoordinate), - utils.getPhiBinIndex(phiCoordinate))} + utils.getPhiBinIndex(phi))} //, montecarloId{ other.montecarloId } { // Nothing to do @@ -56,12 +56,12 @@ Cluster::Cluster(const int layerIndex, const float3& primaryVertex, const IndexT : xCoordinate{other.xCoordinate}, yCoordinate{other.yCoordinate}, zCoordinate{other.zCoordinate}, - phiCoordinate{getNormalizedPhiCoordinate( - calculatePhiCoordinate(xCoordinate - primaryVertex.x, yCoordinate - primaryVertex.y))}, - rCoordinate{calculateRCoordinate(xCoordinate - primaryVertex.x, yCoordinate - primaryVertex.y)}, + phi{getNormalizedPhi( + computePhi(xCoordinate - primaryVertex.x, yCoordinate - primaryVertex.y))}, + radius{hypot(xCoordinate - primaryVertex.x, yCoordinate - primaryVertex.y)}, clusterId{other.clusterId}, indexTableBinIndex{utils.getBinIndex(utils.getZBinIndex(layerIndex, zCoordinate), - utils.getPhiBinIndex(phiCoordinate))} + utils.getPhiBinIndex(phi))} { // Nothing to do } @@ -71,12 +71,12 @@ void Cluster::Init(const int layerIndex, const float3& primaryVertex, const Inde xCoordinate = other.xCoordinate; yCoordinate = other.yCoordinate; zCoordinate = other.zCoordinate; - phiCoordinate = getNormalizedPhiCoordinate( - calculatePhiCoordinate(xCoordinate - primaryVertex.x, yCoordinate - primaryVertex.y)); - rCoordinate = calculateRCoordinate(xCoordinate - primaryVertex.x, yCoordinate - primaryVertex.y); + phi = getNormalizedPhi( + computePhi(xCoordinate - primaryVertex.x, yCoordinate - primaryVertex.y)); + radius = hypot(xCoordinate - primaryVertex.x, yCoordinate - primaryVertex.y); clusterId = other.clusterId; indexTableBinIndex = utils.getBinIndex(utils.getZBinIndex(layerIndex, zCoordinate), - utils.getPhiBinIndex(phiCoordinate)); + utils.getPhiBinIndex(phi)); } TrackingFrameInfo::TrackingFrameInfo(float x, float y, float z, float xTF, float alpha, GPUArray&& posTF, diff --git a/Detectors/ITSMFT/ITS/tracking/src/IOUtils.cxx b/Detectors/ITSMFT/ITS/tracking/src/IOUtils.cxx index fd18501491562..b03374b2c5ee4 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/IOUtils.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/IOUtils.cxx @@ -265,7 +265,7 @@ std::vector> ioutils::loadLabels(const int events std::ifstream inputStream{}; std::string line{}; int monteCarloId{}, pdgCode{}, numberOfClusters{}; - float transverseMomentum{}, phiCoordinate{}, pseudorapidity{}; + float transverseMomentum{}, phi{}, pseudorapidity{}; labelsMap.reserve(eventsNum); @@ -285,12 +285,12 @@ std::vector> ioutils::loadLabels(const int events } else { - if (inputStringStream >> transverseMomentum >> phiCoordinate >> pseudorapidity >> pdgCode >> numberOfClusters) { + if (inputStringStream >> transverseMomentum >> phi >> pseudorapidity >> pdgCode >> numberOfClusters) { if (std::abs(pdgCode) == constants::pdgcodes::PionCode && numberOfClusters == 7) { currentEventLabelsMap.emplace(std::piecewise_construct, std::forward_as_tuple(monteCarloId), - std::forward_as_tuple(monteCarloId, transverseMomentum, phiCoordinate, + std::forward_as_tuple(monteCarloId, transverseMomentum, phi, pseudorapidity, pdgCode, numberOfClusters)); } } diff --git a/Detectors/ITSMFT/ITS/tracking/src/Label.cxx b/Detectors/ITSMFT/ITS/tracking/src/Label.cxx index 16812e4c97ac1..e195318828f51 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Label.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Label.cxx @@ -23,7 +23,7 @@ namespace its Label::Label(const int mcId, const float pT, const float phi, const float eta, const int pdg, const int ncl) : monteCarloId{mcId}, transverseMomentum{pT}, - phiCoordinate{phi}, + phi{phi}, pseudorapidity{eta}, pdgCode{pdg}, numberOfClusters{ncl} @@ -33,7 +33,7 @@ Label::Label(const int mcId, const float pT, const float phi, const float eta, c std::ostream& operator<<(std::ostream& outputStream, const Label& label) { - outputStream << label.monteCarloId << "\t" << label.transverseMomentum << "\t" << label.phiCoordinate << "\t" + outputStream << label.monteCarloId << "\t" << label.transverseMomentum << "\t" << label.phi << "\t" << label.pseudorapidity << "\t" << label.pdgCode << "\t" << label.numberOfClusters; return outputStream; diff --git a/Detectors/ITSMFT/ITS/tracking/src/PrimaryVertexContext.cxx b/Detectors/ITSMFT/ITS/tracking/src/PrimaryVertexContext.cxx index 0b70063e550cb..22c259d7390fd 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/PrimaryVertexContext.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/PrimaryVertexContext.cxx @@ -75,12 +75,12 @@ void PrimaryVertexContext::initialise(const MemoryParameters& memParam, const Tr ClusterHelper& h = cHelper[iCluster]; float x = c.xCoordinate - mPrimaryVertex.x; float y = c.yCoordinate - mPrimaryVertex.y; - float phi = math_utils::calculatePhiCoordinate(x, y); + float phi = math_utils::computePhi(x, y); const int zBin{mIndexTableUtils.getZBinIndex(iLayer, c.zCoordinate)}; int bin = mIndexTableUtils.getBinIndex(zBin, mIndexTableUtils.getPhiBinIndex(phi)); CA_DEBUGGER(assert(zBin > 0)); h.phi = phi; - h.r = math_utils::calculateRCoordinate(x, y); + h.r = std::hypot(x, y); mMinR[iLayer] = o2::gpu::GPUCommonMath::Min(h.r, mMinR[iLayer]); mMaxR[iLayer] = o2::gpu::GPUCommonMath::Max(h.r, mMaxR[iLayer]); h.bin = bin; @@ -97,8 +97,8 @@ void PrimaryVertexContext::initialise(const MemoryParameters& memParam, const Tr ClusterHelper& h = cHelper[iCluster]; Cluster& c = mClusters[iLayer][lutPerBin[h.bin] + h.ind]; c = currentLayer[iCluster]; - c.phiCoordinate = h.phi; - c.rCoordinate = h.r; + c.phi = h.phi; + c.radius = h.r; c.indexTableBinIndex = h.bin; } diff --git a/Detectors/ITSMFT/ITS/tracking/src/StandaloneDebugger.cxx b/Detectors/ITSMFT/ITS/tracking/src/StandaloneDebugger.cxx index 7be852282841f..662fc1a0f6ff4 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/StandaloneDebugger.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/StandaloneDebugger.cxx @@ -63,7 +63,7 @@ void StandaloneDebugger::fillCombinatoricsTree(std::array, (*mTreeStream) << "combinatorics01" << "tanLambda=" << combination.tanLambda - << "phi=" << combination.phiCoordinate + << "phi=" << combination.phi << "c0z=" << c0z << "c1z=" << c1z << "isValidated=" << isValidated @@ -81,7 +81,7 @@ void StandaloneDebugger::fillCombinatoricsTree(std::array, (*mTreeStream) << "combinatorics12" << "tanLambda=" << combination.tanLambda - << "phi=" << combination.phiCoordinate + << "phi=" << combination.phi << "c1z=" << c1z << "c2z=" << c2z << "isValidated=" << isValidated @@ -105,7 +105,7 @@ void StandaloneDebugger::fillTrackletSelectionTree(std::arraygetClusterLabels(1, clusters[1][comb01[trackletPair[0]].secondClusterIndex].clusterId); o2::MCCompLabel lblClus2 = event->getClusterLabels(2, clusters[2][comb12[trackletPair[1]].secondClusterIndex].clusterId); unsigned char isValidated{lblClus0.compare(lblClus1) == 1 && lblClus0.compare(lblClus2) == 1}; - float deltaPhi{comb01[trackletPair[0]].phiCoordinate - comb12[trackletPair[1]].phiCoordinate}; + float deltaPhi{comb01[trackletPair[0]].phi - comb12[trackletPair[1]].phi}; float deltaTanLambda{comb01[trackletPair[0]].tanLambda - comb12[trackletPair[1]].tanLambda}; mTreeStream->GetDirectory()->cd(); // in case of existing other open files (*mTreeStream) @@ -115,13 +115,13 @@ void StandaloneDebugger::fillTrackletSelectionTree(std::array + +namespace +{ +struct ClusterHelper { + float phi; + float r; + int bin; + int ind; +}; +} // namespace + +namespace o2 +{ +namespace its +{ + +constexpr float DefClusErrorRow = o2::itsmft::SegmentationAlpide::PitchRow * 0.5; +constexpr float DefClusErrorCol = o2::itsmft::SegmentationAlpide::PitchCol * 0.5; +constexpr float DefClusError2Row = DefClusErrorRow * DefClusErrorRow; +constexpr float DefClusError2Col = DefClusErrorCol * DefClusErrorCol; + +TimeFrame::TimeFrame(int nLayers) +{ + mMinR.resize(nLayers, 10000.); + mMaxR.resize(nLayers, -1.); + mClusters.resize(nLayers); + mUnsortedClusters.resize(nLayers); + mTrackingFrameInfo.resize(nLayers); + mClusterLabels.resize(nLayers); + mClusterExternalIndices.resize(nLayers); + mUsedClusters.resize(nLayers); + mROframesClusters.resize(nLayers, {0}); ///TBC: if resetting the timeframe is required, then this has to be done +} + +void TimeFrame::addPrimaryVertices(const std::vector& vertices) +{ + for (const auto& vertex : vertices) { + mPrimaryVertices.emplace_back(vertex); + const int w{vertex.getNContributors()}; + mBeamPos[0] = (mBeamPos[0] * mBeamPosWeight + vertex.getX() * w) / (mBeamPosWeight + w); + mBeamPos[1] = (mBeamPos[1] * mBeamPosWeight + vertex.getY() * w) / (mBeamPosWeight + w); + mBeamPosWeight += w; + } + mROframesPV.push_back(mPrimaryVertices.size()); +} + +int TimeFrame::loadROFrameData(const o2::itsmft::ROFRecord& rof, gsl::span clusters, + const dataformats::MCTruthContainer* mcLabels) +{ + GeometryTGeo* geom = GeometryTGeo::Instance(); + geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G)); + int clusterId{0}; + + auto first = rof.getFirstEntry(); + auto clusters_in_frame = rof.getROFData(clusters); + for (auto& c : clusters_in_frame) { + int layer = geom->getLayer(c.getSensorID()); + + /// Clusters are stored in the tracking frame + auto xyz = c.getXYZGloRot(*geom); + addTrackingFrameInfoToLayer(layer, xyz.x(), xyz.y(), xyz.z(), c.getX(), geom->getSensorRefAlpha(c.getSensorID()), + std::array{c.getY(), c.getZ()}, + std::array{c.getSigmaY2(), c.getSigmaYZ(), c.getSigmaZ2()}); + + /// Rotate to the global frame + addClusterToLayer(layer, xyz.x(), xyz.y(), xyz.z(), mUnsortedClusters[layer].size()); + if (mcLabels) { + addClusterLabelToLayer(layer, *(mcLabels->getLabels(first + clusterId).begin())); + } + addClusterExternalIndexToLayer(layer, first + clusterId); + clusterId++; + } + + for (unsigned int iL{0}; iL < mUnsortedClusters.size(); ++iL) { + mROframesClusters[iL].push_back(mUnsortedClusters[iL].size()); + } + mNrof++; + return clusters_in_frame.size(); +} + +int TimeFrame::loadROFrameData(gsl::span rofs, gsl::span clusters, gsl::span::iterator& pattIt, const itsmft::TopologyDictionary& dict, const dataformats::MCTruthContainer* mcLabels) +{ + GeometryTGeo* geom = GeometryTGeo::Instance(); + geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G)); + + mNrof = 0; + for (int clusterId{0}; clusterId < clusters.size() && mNrof < rofs.size(); ++clusterId) { + auto& c = clusters[clusterId]; + + int layer = geom->getLayer(c.getSensorID()); + + auto pattID = c.getPatternID(); + o2::math_utils::Point3D locXYZ; + float sigmaY2 = DefClusError2Row, sigmaZ2 = DefClusError2Col, sigmaYZ = 0; //Dummy COG errors (about half pixel size) + if (pattID != itsmft::CompCluster::InvalidPatternID) { + sigmaY2 = dict.getErr2X(pattID); + sigmaZ2 = dict.getErr2Z(pattID); + if (!dict.isGroup(pattID)) { + locXYZ = dict.getClusterCoordinates(c); + } else { + o2::itsmft::ClusterPattern patt(pattIt); + locXYZ = dict.getClusterCoordinates(c, patt); + } + } else { + o2::itsmft::ClusterPattern patt(pattIt); + locXYZ = dict.getClusterCoordinates(c, patt, false); + } + auto sensorID = c.getSensorID(); + // Inverse transformation to the local --> tracking + auto trkXYZ = geom->getMatrixT2L(sensorID) ^ locXYZ; + // Transformation to the local --> global + auto gloXYZ = geom->getMatrixL2G(sensorID) * locXYZ; + + addTrackingFrameInfoToLayer(layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), trkXYZ.x(), geom->getSensorRefAlpha(sensorID), + std::array{trkXYZ.y(), trkXYZ.z()}, + std::array{sigmaY2, sigmaYZ, sigmaZ2}); + + /// Rotate to the global frame + addClusterToLayer(layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), mUnsortedClusters[layer].size()); + if (mcLabels) { + addClusterLabelToLayer(layer, *(mcLabels->getLabels(clusterId).begin())); + } + addClusterExternalIndexToLayer(layer, clusterId); + + if (clusterId == rofs[mNrof].getFirstEntry() + rofs[mNrof].getNEntries() - 1) { + for (unsigned int iL{0}; iL < mUnsortedClusters.size(); ++iL) { + mROframesClusters[iL].push_back(mUnsortedClusters[iL].size()); + } + mNrof++; + } + } + return mNrof; +} + +int TimeFrame::getTotalClusters() const +{ + size_t totalClusters{0}; + for (auto& clusters : mUnsortedClusters) { + totalClusters += clusters.size(); + } + return int(totalClusters); +} + +void TimeFrame::initialise(const int iteration, const MemoryParameters& memParam, const TrackingParameters& trkParam) +{ + if (iteration == 0) { + mTracks.clear(); + mTracksLabel.clear(); + mTracks.resize(mNrof); + mTracksLabel.resize(mNrof); + mCells.resize(trkParam.CellsPerRoad()); + mCellsLookupTable.resize(trkParam.CellsPerRoad() - 1); + mCellsNeighbours.resize(trkParam.CellsPerRoad() - 1); + mCellLabels.resize(trkParam.CellsPerRoad()); + mTracklets.resize(trkParam.TrackletsPerRoad()); + mTrackletLabels.resize(trkParam.TrackletsPerRoad()); + mTrackletsLookupTable.resize(trkParam.CellsPerRoad()); + mIndexTableUtils.setTrackingParameters(trkParam); + + for (unsigned int iLayer{0}; iLayer < mClusters.size(); ++iLayer) { + if (mClusters[iLayer].size()) { + continue; + } + mClusters[iLayer].clear(); + mClusters[iLayer].resize(mUnsortedClusters[iLayer].size()); + mUsedClusters[iLayer].clear(); + mUsedClusters[iLayer].resize(mUnsortedClusters[iLayer].size(), false); + } + + mIndexTables.resize(mNrof); + std::vector cHelper; + for (int rof{0}; rof < mNrof; ++rof) { + mIndexTables[rof].resize(trkParam.TrackletsPerRoad(), std::vector(trkParam.ZBins * trkParam.PhiBins + 1, 0)); + for (int iLayer{0}; iLayer < trkParam.NLayers; ++iLayer) { + std::vector clsPerBin(trkParam.PhiBins * trkParam.ZBins, 0); + + const auto unsortedClusters{getUnsortedClustersOnLayer(rof, iLayer)}; + const int clustersNum{static_cast(unsortedClusters.size())}; + + cHelper.clear(); + cHelper.resize(clustersNum); + + for (int iCluster{0}; iCluster < clustersNum; ++iCluster) { + const Cluster& c = unsortedClusters[iCluster]; + ClusterHelper& h = cHelper[iCluster]; + float x = c.xCoordinate - mBeamPos[0]; + float y = c.yCoordinate - mBeamPos[1]; + float phi = math_utils::computePhi(x, y); + const int zBin{mIndexTableUtils.getZBinIndex(iLayer, c.zCoordinate)}; + int bin = mIndexTableUtils.getBinIndex(zBin, mIndexTableUtils.getPhiBinIndex(phi)); + h.phi = phi; + h.r = math_utils::hypot(x, y); + mMinR[iLayer] = o2::gpu::GPUCommonMath::Min(h.r, mMinR[iLayer]); + mMaxR[iLayer] = o2::gpu::GPUCommonMath::Max(h.r, mMaxR[iLayer]); + h.bin = bin; + h.ind = clsPerBin[bin]++; + } + std::vector lutPerBin(clsPerBin.size()); + lutPerBin[0] = 0; + for (unsigned int iB{1}; iB < lutPerBin.size(); ++iB) { + lutPerBin[iB] = lutPerBin[iB - 1] + clsPerBin[iB - 1]; + } + + auto clusters2beSorted{getClustersOnLayer(rof, iLayer)}; + for (int iCluster{0}; iCluster < clustersNum; ++iCluster) { + const ClusterHelper& h = cHelper[iCluster]; + + Cluster& c = clusters2beSorted[lutPerBin[h.bin] + h.ind]; + c = unsortedClusters[iCluster]; + c.phi = h.phi; + c.radius = h.r; + c.indexTableBinIndex = h.bin; + } + + if (iLayer > 0) { + for (unsigned int iB{0}; iB < clsPerBin.size(); ++iB) { + mIndexTables[rof][iLayer - 1][iB] = lutPerBin[iB]; + } + for (auto iB{clsPerBin.size()}; iB < (int)mIndexTables[rof][iLayer - 1].size(); iB++) { + mIndexTables[rof][iLayer - 1][iB] = clustersNum; + } + } + } + } + } + + mRoads.clear(); + + for (unsigned int iLayer{0}; iLayer < mTracklets.size(); ++iLayer) { + mTracklets[iLayer].clear(); + mTrackletLabels[iLayer].clear(); + if (iLayer < mCells.size()) { + mCells[iLayer].clear(); + mTrackletsLookupTable[iLayer].clear(); + mTrackletsLookupTable[iLayer].resize(mClusters[iLayer + 1].size(), 0); + mCellLabels[iLayer].clear(); + } + + if (iLayer < mCells.size() - 1) { + mCellsLookupTable[iLayer].clear(); + mCellsNeighbours[iLayer].clear(); + } + } +} + +void TimeFrame::checkTrackletLUTs() +{ + for (uint32_t iLayer{0}; iLayer < getTracklets().size(); ++iLayer) { + int prev{-1}; + int count{0}; + for (uint32_t iTracklet{0}; iTracklet < getTracklets()[iLayer].size(); ++iTracklet) { + auto& trk = getTracklets()[iLayer][iTracklet]; + int currentId{trk.firstClusterIndex}; + if (currentId < prev) { + std::cout << "First Cluster Index not increasing monotonically on L:T:ID:Prev " << iLayer << "\t" << iTracklet << "\t" << currentId << "\t" << prev << std::endl; + } else if (currentId == prev) { + count++; + } else { + if (iLayer > 0) { + auto& lut{getTrackletsLookupTable()[iLayer - 1]}; + if (count != lut[prev + 1] - lut[prev]) { + std::cout << "LUT count broken " << iLayer - 1 << "\t" << prev << "\t" << count << "\t" << lut[prev + 1] << "\t" << lut[prev] << std::endl; + } + } + count = 1; + } + prev = currentId; + if (iLayer > 0) { + auto& lut{getTrackletsLookupTable()[iLayer - 1]}; + if (iTracklet >= lut[currentId + 1] || iTracklet < lut[currentId]) { + std::cout << "LUT broken: " << iLayer - 1 << "\t" << currentId << "\t" << iTracklet << std::endl; + } + } + } + } +} + +void TimeFrame::printTrackletLUTonLayer(int i) +{ + std::cout << "--------" << std::endl + << "Tracklet LUT " << i << std::endl; + for (int j : mTrackletsLookupTable[i]) { + std::cout << j << "\t"; + } + std::cout << "\n--------" << std::endl + << std::endl; +} + +void TimeFrame::printCellLUTonLayer(int i) +{ + std::cout << "--------" << std::endl + << "Cell LUT " << i << std::endl; + for (int j : mCellsLookupTable[i]) { + std::cout << j << "\t"; + } + std::cout << "\n--------" << std::endl + << std::endl; +} + +void TimeFrame::printTrackletLUTs() +{ + for (unsigned int i{0}; i < mTrackletsLookupTable.size(); ++i) { + printTrackletLUTonLayer(i); + } +} + +void TimeFrame::printCellLUTs() +{ + for (unsigned int i{0}; i < mCellsLookupTable.size(); ++i) { + printCellLUTonLayer(i); + } +} + +void TimeFrame::printVertices() +{ + std::cout << "Vertices in ROF (nROF = " << mNrof << ", lut size = " << mROframesPV.size() << ")" << std::endl; + for (unsigned int iR{0}; iR < mROframesPV.size(); ++iR) { + std::cout << mROframesPV[iR] << "\t"; + } + std::cout << "\n\n Vertices:" << std::endl; + for (unsigned int iV{0}; iV < mPrimaryVertices.size(); ++iV) { + std::cout << mPrimaryVertices[iV].getX() << "\t" << mPrimaryVertices[iV].getY() << "\t" << mPrimaryVertices[iV].getZ() << std::endl; + } + std::cout << "--------" << std::endl; +} + +void TimeFrame::printROFoffsets() +{ + std::cout << "--------" << std::endl; + for (unsigned int iLayer{0}; iLayer < mROframesClusters.size(); ++iLayer) { + std::cout << "Layer " << iLayer << std::endl; + for (auto value : mROframesClusters[iLayer]) { + std::cout << value << "\t"; + } + std::cout << std::endl; + } +} + +} // namespace its +} // namespace o2 diff --git a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx index 2eb5dfbb9adfd..a45b16ccf9670 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx @@ -29,6 +29,7 @@ #include #include #include +#include namespace o2 { @@ -41,7 +42,7 @@ Tracker::Tracker(o2::its::TrackerTraits* traits) mTrkParams.resize(1); mMemParams.resize(1); mTraits = traits; - mPrimaryVertexContext = mTraits->getPrimaryVertexContext(); + mTimeFrame = mTraits->getTimeFrame(); #ifdef CA_DEBUG mDebugger = new StandaloneDebugger("dbg_ITSTrackerCPU.root"); #endif @@ -55,48 +56,32 @@ Tracker::~Tracker() Tracker::~Tracker() = default; #endif -void Tracker::clustersToTracks(const ROframe& event, std::ostream& timeBenchmarkOutputStream) +void Tracker::clustersToTracks(std::function logger) { - const int verticesNum = event.getPrimaryVerticesNum(); - mTracks.clear(); - mTrackLabels.clear(); - - for (int iVertex = 0; iVertex < verticesNum; ++iVertex) { - - float total{0.f}; - - for (int iteration = 0; iteration < mTrkParams.size(); ++iteration) { - - int numCls = 0; - for (unsigned int iLayer{0}; iLayer < mTrkParams[iteration].NLayers; ++iLayer) { - numCls += event.getClusters()[iLayer].size(); - } - if (numCls < mTrkParams[iteration].MinTrackLength) { - continue; - } + double total{0}; + for (int iteration = 0; iteration < mTrkParams.size(); ++iteration) { + mTraits->UpdateTrackingParameters(mTrkParams[iteration]); + + total += evaluateTask(&Tracker::initialiseTimeFrame, "Context initialisation", + logger, iteration, mMemParams[iteration], mTrkParams[iteration]); + total += evaluateTask(&Tracker::computeTracklets, "Tracklet finding", logger); + total += evaluateTask(&Tracker::computeCells, "Cell finding", logger); + total += evaluateTask(&Tracker::findCellsNeighbours, "Neighbour finding", logger, iteration); + total += evaluateTask(&Tracker::findRoads, "Road finding", logger, iteration); + total += evaluateTask(&Tracker::findTracks, "Track finding", logger); + } - mTraits->UpdateTrackingParameters(mTrkParams[iteration]); - /// Ugly hack -> Unifiy float3 definition in CPU and CUDA/HIP code - int pass = iteration + iVertex; /// Do not reinitialise the context if we analyse pile-up events - std::array pV = {event.getPrimaryVertex(iVertex).x, event.getPrimaryVertex(iVertex).y, event.getPrimaryVertex(iVertex).z}; - total += evaluateTask(&Tracker::initialisePrimaryVertexContext, "Context initialisation", - timeBenchmarkOutputStream, mMemParams[iteration], mTrkParams[iteration], event.getClusters(), pV, pass); - total += evaluateTask(&Tracker::computeTracklets, "Tracklet finding", timeBenchmarkOutputStream); - total += evaluateTask(&Tracker::computeCells, "Cell finding", timeBenchmarkOutputStream); - total += evaluateTask(&Tracker::findCellsNeighbours, "Neighbour finding", timeBenchmarkOutputStream, iteration); - total += evaluateTask(&Tracker::findRoads, "Road finding", timeBenchmarkOutputStream, iteration); - total += evaluateTask(&Tracker::findTracks, "Track finding", timeBenchmarkOutputStream, event); - } - if (constants::DoTimeBenchmarks && fair::Logger::Logging(fair::Severity::info)) { - timeBenchmarkOutputStream << std::setw(2) << " - " - << "Vertex processing completed in: " << total << "ms" << std::endl; - } + std::stringstream sstream; + if (constants::DoTimeBenchmarks) { + sstream << std::setw(2) << " - " + << "Vertex processing completed in: " << total << "ms" << std::endl; } - if (event.hasMCinformation()) { - computeTracksMClabels(event); - } else { - rectifyClusterIndices(event); + logger(sstream.str()); + + if (mTimeFrame->hasMCinformation()) { + computeTracksMClabels(); } + rectifyClusterIndices(); } void Tracker::computeTracklets() @@ -113,53 +98,48 @@ void Tracker::findCellsNeighbours(int& iteration) { for (int iLayer{0}; iLayer < mTrkParams[iteration].CellsPerRoad() - 1; ++iLayer) { - if (mPrimaryVertexContext->getCells()[iLayer + 1].empty() || - mPrimaryVertexContext->getCellsLookupTable()[iLayer].empty()) { + if (mTimeFrame->getCells()[iLayer + 1].empty() || + mTimeFrame->getCellsLookupTable()[iLayer].empty()) { continue; } - int layerCellsNum{static_cast(mPrimaryVertexContext->getCells()[iLayer].size())}; - const int nextLayerCellsNum{static_cast(mPrimaryVertexContext->getCells()[iLayer + 1].size())}; - mPrimaryVertexContext->getCellsNeighbours()[iLayer].resize(nextLayerCellsNum); + int layerCellsNum{static_cast(mTimeFrame->getCells()[iLayer].size())}; + const int nextLayerCellsNum{static_cast(mTimeFrame->getCells()[iLayer + 1].size())}; + mTimeFrame->getCellsNeighbours()[iLayer].resize(nextLayerCellsNum); for (int iCell{0}; iCell < layerCellsNum; ++iCell) { - const Cell& currentCell{mPrimaryVertexContext->getCells()[iLayer][iCell]}; + const Cell& currentCell{mTimeFrame->getCells()[iLayer][iCell]}; const int nextLayerTrackletIndex{currentCell.getSecondTrackletIndex()}; - const int nextLayerFirstCellIndex{mPrimaryVertexContext->getCellsLookupTable()[iLayer][nextLayerTrackletIndex]}; - if (nextLayerFirstCellIndex != constants::its::UnusedIndex && - mPrimaryVertexContext->getCells()[iLayer + 1][nextLayerFirstCellIndex].getFirstTrackletIndex() == - nextLayerTrackletIndex) { - - for (int iNextLayerCell{nextLayerFirstCellIndex}; iNextLayerCell < nextLayerCellsNum; ++iNextLayerCell) { - - Cell& nextCell{mPrimaryVertexContext->getCells()[iLayer + 1][iNextLayerCell]}; - if (nextCell.getFirstTrackletIndex() != nextLayerTrackletIndex) { - break; - } + const int nextLayerFirstCellIndex{mTimeFrame->getCellsLookupTable()[iLayer][nextLayerTrackletIndex]}; + const int nextLayerLastCellIndex{mTimeFrame->getCellsLookupTable()[iLayer][nextLayerTrackletIndex + 1]}; + for (int iNextLayerCell{nextLayerFirstCellIndex}; iNextLayerCell < nextLayerLastCellIndex; ++iNextLayerCell) { - const float3 currentCellNormalVector{currentCell.getNormalVectorCoordinates()}; - const float3 nextCellNormalVector{nextCell.getNormalVectorCoordinates()}; - const float3 normalVectorsDeltaVector{currentCellNormalVector.x - nextCellNormalVector.x, - currentCellNormalVector.y - nextCellNormalVector.y, - currentCellNormalVector.z - nextCellNormalVector.z}; + Cell& nextCell{mTimeFrame->getCells()[iLayer + 1][iNextLayerCell]}; + if (nextCell.getFirstTrackletIndex() != nextLayerTrackletIndex) { + break; + } - const float deltaNormalVectorsModulus{(normalVectorsDeltaVector.x * normalVectorsDeltaVector.x) + - (normalVectorsDeltaVector.y * normalVectorsDeltaVector.y) + - (normalVectorsDeltaVector.z * normalVectorsDeltaVector.z)}; - const float deltaCurvature{std::abs(currentCell.getCurvature() - nextCell.getCurvature())}; + const float3 currentCellNormalVector{currentCell.getNormalVectorCoordinates()}; + const float3 nextCellNormalVector{nextCell.getNormalVectorCoordinates()}; + const float3 normalVectorsDeltaVector{currentCellNormalVector.x - nextCellNormalVector.x, + currentCellNormalVector.y - nextCellNormalVector.y, + currentCellNormalVector.z - nextCellNormalVector.z}; - if (deltaNormalVectorsModulus < mTrkParams[iteration].NeighbourMaxDeltaN[iLayer] && - deltaCurvature < mTrkParams[iteration].NeighbourMaxDeltaCurvature[iLayer]) { + const float deltaNormalVectorsModulus{(normalVectorsDeltaVector.x * normalVectorsDeltaVector.x) + + (normalVectorsDeltaVector.y * normalVectorsDeltaVector.y) + + (normalVectorsDeltaVector.z * normalVectorsDeltaVector.z)}; + const float deltaCurvature{std::abs(currentCell.getCurvature() - nextCell.getCurvature())}; - mPrimaryVertexContext->getCellsNeighbours()[iLayer][iNextLayerCell].push_back(iCell); + if (deltaNormalVectorsModulus < mTrkParams[iteration].NeighbourMaxDeltaN[iLayer] && + deltaCurvature < mTrkParams[iteration].NeighbourMaxDeltaCurvature[iLayer]) { - const int currentCellLevel{currentCell.getLevel()}; + mTimeFrame->getCellsNeighbours()[iLayer][iNextLayerCell].push_back(iCell); - if (currentCellLevel >= nextCell.getLevel()) { + const int currentCellLevel{currentCell.getLevel()}; - nextCell.setLevel(currentCellLevel + 1); - } + if (currentCellLevel >= nextCell.getLevel()) { + nextCell.setLevel(currentCellLevel + 1); } } } @@ -170,22 +150,22 @@ void Tracker::findCellsNeighbours(int& iteration) void Tracker::findRoads(int& iteration) { for (int iLevel{mTrkParams[iteration].CellsPerRoad()}; iLevel >= mTrkParams[iteration].CellMinimumLevel(); --iLevel) { - CA_DEBUGGER(int nRoads = -mPrimaryVertexContext->getRoads().size()); + CA_DEBUGGER(int nRoads = -mTimeFrame->getRoads().size()); const int minimumLevel{iLevel - 1}; for (int iLayer{mTrkParams[iteration].CellsPerRoad() - 1}; iLayer >= minimumLevel; --iLayer) { - const int levelCellsNum{static_cast(mPrimaryVertexContext->getCells()[iLayer].size())}; + const int levelCellsNum{static_cast(mTimeFrame->getCells()[iLayer].size())}; for (int iCell{0}; iCell < levelCellsNum; ++iCell) { - Cell& currentCell{mPrimaryVertexContext->getCells()[iLayer][iCell]}; + Cell& currentCell{mTimeFrame->getCells()[iLayer][iCell]}; if (currentCell.getLevel() != iLevel) { continue; } - mPrimaryVertexContext->getRoads().emplace_back(iLayer, iCell); + mTimeFrame->getRoads().emplace_back(iLayer, iCell); /// For 3 clusters roads (useful for cascades and hypertriton) we just store the single cell /// and we do not do the candidate tree traversal @@ -194,13 +174,13 @@ void Tracker::findRoads(int& iteration) } const int cellNeighboursNum{static_cast( - mPrimaryVertexContext->getCellsNeighbours()[iLayer - 1][iCell].size())}; + mTimeFrame->getCellsNeighbours()[iLayer - 1][iCell].size())}; bool isFirstValidNeighbour = true; for (int iNeighbourCell{0}; iNeighbourCell < cellNeighboursNum; ++iNeighbourCell) { - const int neighbourCellId = mPrimaryVertexContext->getCellsNeighbours()[iLayer - 1][iCell][iNeighbourCell]; - const Cell& neighbourCell = mPrimaryVertexContext->getCells()[iLayer - 1][neighbourCellId]; + const int neighbourCellId = mTimeFrame->getCellsNeighbours()[iLayer - 1][iCell][iNeighbourCell]; + const Cell& neighbourCell = mTimeFrame->getCells()[iLayer - 1][neighbourCellId]; if (iLevel - 1 != neighbourCell.getLevel()) { continue; @@ -212,7 +192,7 @@ void Tracker::findRoads(int& iteration) } else { - mPrimaryVertexContext->getRoads().emplace_back(iLayer, iCell); + mTimeFrame->getRoads().emplace_back(iLayer, iCell); } traverseCellsTree(neighbourCellId, iLayer - 1); @@ -223,39 +203,37 @@ void Tracker::findRoads(int& iteration) } } #ifdef CA_DEBUG - nRoads += mPrimaryVertexContext->getRoads().size(); - std::cout << "+++ Roads with " << iLevel + 2 << " clusters: " << nRoads << " / " << mPrimaryVertexContext->getRoads().size() << std::endl; + nRoads += mTimeFrame->getRoads().size(); + std::cout << "+++ Roads with " << iLevel + 2 << " clusters: " << nRoads << " / " << mTimeFrame->getRoads().size() << std::endl; #endif } } -void Tracker::findTracks(const ROframe& event) +void Tracker::findTracks() { - mTracks.reserve(mTracks.capacity() + mPrimaryVertexContext->getRoads().size()); std::vector tracks; - tracks.reserve(mPrimaryVertexContext->getRoads().size()); + tracks.reserve(mTimeFrame->getRoads().size()); -#ifdef CA_DEBUG - std::vector roadCounters(mTrkParams[0].NLayers - 3, 0); - std::vector fitCounters(mTrkParams[0].NLayers - 3, 0); - std::vector backpropagatedCounters(mTrkParams[0].NLayers - 3, 0); - std::vector refitCounters(mTrkParams[0].NLayers - 3, 0); - std::vector nonsharingCounters(mTrkParams[0].NLayers - 3, 0); -#endif - - for (auto& road : mPrimaryVertexContext->getRoads()) { + for (auto& road : mTimeFrame->getRoads()) { std::vector clusters(mTrkParams[0].NLayers, constants::its::UnusedIndex); int lastCellLevel = constants::its::UnusedIndex; CA_DEBUGGER(int nClusters = 2); + int firstTracklet{constants::its::UnusedIndex}; + std::vector tracklets(mTrkParams[0].TrackletsPerRoad(), constants::its::UnusedIndex); for (int iCell{0}; iCell < mTrkParams[0].CellsPerRoad(); ++iCell) { const int cellIndex = road[iCell]; if (cellIndex == constants::its::UnusedIndex) { continue; } else { - clusters[iCell] = mPrimaryVertexContext->getCells()[iCell][cellIndex].getFirstClusterIndex(); - clusters[iCell + 1] = mPrimaryVertexContext->getCells()[iCell][cellIndex].getSecondClusterIndex(); - clusters[iCell + 2] = mPrimaryVertexContext->getCells()[iCell][cellIndex].getThirdClusterIndex(); + if (firstTracklet == constants::its::UnusedIndex) { + firstTracklet = iCell; + } + tracklets[iCell] = mTimeFrame->getCells()[iCell][cellIndex].getFirstTrackletIndex(); + tracklets[iCell + 1] = mTimeFrame->getCells()[iCell][cellIndex].getSecondTrackletIndex(); + clusters[iCell] = mTimeFrame->getCells()[iCell][cellIndex].getFirstClusterIndex(); + clusters[iCell + 1] = mTimeFrame->getCells()[iCell][cellIndex].getSecondClusterIndex(); + clusters[iCell + 2] = mTimeFrame->getCells()[iCell][cellIndex].getThirdClusterIndex(); assert(clusters[iCell] != constants::its::UnusedIndex && clusters[iCell + 1] != constants::its::UnusedIndex && clusters[iCell + 2] != constants::its::UnusedIndex); @@ -264,6 +242,24 @@ void Tracker::findTracks(const ROframe& event) } } + CA_DEBUGGER(assert(nClusters >= mTrkParams[0].MinTrackLength)); + int count{1}; + unsigned short rof{mTimeFrame->getTracklets()[firstTracklet][tracklets[firstTracklet]].rof[0]}; + for (int iT = firstTracklet; iT < 6; ++iT) { + if (tracklets[iT] == constants::its::UnusedIndex) { + continue; + } + if (rof == mTimeFrame->getTracklets()[iT][tracklets[iT]].rof[1]) { + count++; + } else { + if (count == 1) { + rof = mTimeFrame->getTracklets()[iT][tracklets[iT]].rof[1]; + } else { + count--; + } + } + } + CA_DEBUGGER(assert(nClusters >= mTrkParams[0].MinTrackLength)); CA_DEBUGGER(roadCounters[nClusters - 4]++); @@ -274,151 +270,96 @@ void Tracker::findTracks(const ROframe& event) /// From primary vertex context index to event index (== the one used as input of the tracking code) for (int iC{0}; iC < clusters.size(); iC++) { if (clusters[iC] != constants::its::UnusedIndex) { - clusters[iC] = mPrimaryVertexContext->getClusters()[iC][clusters[iC]].clusterId; + clusters[iC] = mTimeFrame->getClusters()[iC][clusters[iC]].clusterId; } } /// Track seed preparation. Clusters are numbered progressively from the outermost to the innermost. - const auto& cluster1_glo = event.getClustersOnLayer(lastCellLevel + 2).at(clusters[lastCellLevel + 2]); - const auto& cluster2_glo = event.getClustersOnLayer(lastCellLevel + 1).at(clusters[lastCellLevel + 1]); - const auto& cluster3_glo = event.getClustersOnLayer(lastCellLevel).at(clusters[lastCellLevel]); + const auto& cluster1_glo = mTimeFrame->getUnsortedClusters()[lastCellLevel + 2].at(clusters[lastCellLevel + 2]); + const auto& cluster2_glo = mTimeFrame->getUnsortedClusters()[lastCellLevel + 1].at(clusters[lastCellLevel + 1]); + const auto& cluster3_glo = mTimeFrame->getUnsortedClusters()[lastCellLevel].at(clusters[lastCellLevel]); - const auto& cluster3_tf = event.getTrackingFrameInfoOnLayer(lastCellLevel).at(clusters[lastCellLevel]); + const auto& cluster3_tf = mTimeFrame->getTrackingFrameInfoOnLayer(lastCellLevel).at(clusters[lastCellLevel]); /// FIXME! TrackITSExt temporaryTrack{buildTrackSeed(cluster1_glo, cluster2_glo, cluster3_glo, cluster3_tf)}; for (size_t iC = 0; iC < clusters.size(); ++iC) { temporaryTrack.setExternalClusterIndex(iC, clusters[iC], clusters[iC] != constants::its::UnusedIndex); } - bool fitSuccess = fitTrack(event, temporaryTrack, mTrkParams[0].NLayers - 4, -1, -1); + bool fitSuccess = fitTrack(temporaryTrack, mTrkParams[0].NLayers - 4, -1, -1); if (!fitSuccess) { continue; } CA_DEBUGGER(fitCounters[nClusters - 4]++); temporaryTrack.resetCovariance(); - fitSuccess = fitTrack(event, temporaryTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].FitIterationMaxChi2[0]); + fitSuccess = fitTrack(temporaryTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].FitIterationMaxChi2[0]); if (!fitSuccess) { continue; } CA_DEBUGGER(backpropagatedCounters[nClusters - 4]++); temporaryTrack.getParamOut() = temporaryTrack; temporaryTrack.resetCovariance(); - fitSuccess = fitTrack(event, temporaryTrack, mTrkParams[0].NLayers - 1, -1, -1, mTrkParams[0].FitIterationMaxChi2[1]); -#ifdef CA_DEBUG - mDebugger->dumpTrackToBranchWithInfo("testBranch", temporaryTrack, event, mPrimaryVertexContext, true); -#endif + fitSuccess = fitTrack(temporaryTrack, mTrkParams[0].NLayers - 1, -1, -1, mTrkParams[0].FitIterationMaxChi2[1], 50.); if (!fitSuccess) { continue; } - CA_DEBUGGER(refitCounters[nClusters - 4]++); + // temporaryTrack.setROFrame(rof); tracks.emplace_back(temporaryTrack); - CA_DEBUGGER(assert(nClusters == temporaryTrack.getNumberOfClusters())); } //mTraits->refitTracks(event.getTrackingFrameInfo(), tracks); std::sort(tracks.begin(), tracks.end(), [](TrackITSExt& track1, TrackITSExt& track2) { return track1.isBetter(track2, 1.e6f); }); -#ifdef CA_DEBUG - // std::array sharingMatrix{0}; - // int prevNclusters = 7; - // auto cumulativeIndex = [](int ncl) -> int { - // constexpr int idx[5] = {0, 5, 11, 18, 26}; - // return idx[ncl - 4]; - // }; - // std::array xcheckCounters{0}; -#endif - for (auto& track : tracks) { - CA_DEBUGGER(int nClusters = 0); int nShared = 0; for (int iLayer{0}; iLayer < mTrkParams[0].NLayers; ++iLayer) { if (track.getClusterIndex(iLayer) == constants::its::UnusedIndex) { continue; } - nShared += int(mPrimaryVertexContext->isClusterUsed(iLayer, track.getClusterIndex(iLayer))); - CA_DEBUGGER(nClusters++); + nShared += int(mTimeFrame->isClusterUsed(iLayer, track.getClusterIndex(iLayer))); } - // #ifdef CA_DEBUG - // assert(nClusters == track.getNumberOfClusters()); - // xcheckCounters[nClusters - 4]++; - // assert(nShared <= nClusters); - // sharingMatrix[cumulativeIndex(nClusters) + nShared]++; - // #endif - if (nShared > mTrkParams[0].ClusterSharing) { continue; } - // #ifdef CA_DEBUG - // nonsharingCounters[nClusters - 4]++; - // assert(nClusters <= prevNclusters); - // prevNclusters = nClusters; - // #endif - + std::array rofs{INT_MAX, INT_MAX, INT_MAX}; for (int iLayer{0}; iLayer < mTrkParams[0].NLayers; ++iLayer) { if (track.getClusterIndex(iLayer) == constants::its::UnusedIndex) { continue; } - mPrimaryVertexContext->markUsedCluster(iLayer, track.getClusterIndex(iLayer)); + mTimeFrame->markUsedCluster(iLayer, track.getClusterIndex(iLayer)); + int currentROF = mTimeFrame->getClusterROF(iLayer, track.getClusterIndex(iLayer)); + for (int iR{0}; iR < 3; ++iR) { + if (rofs[iR] == INT_MAX) { + rofs[iR] = currentROF; + } + if (rofs[iR] == currentROF) { + break; + } + } + } + if (rofs[2] != INT_MAX) { + continue; + } + if (rofs[1] != INT_MAX) { + track.setNextROFbit(); } - mTracks.emplace_back(track); + mTimeFrame->getTracks(std::min(rofs[0], rofs[1])).emplace_back(track); } - -#ifdef CA_DEBUG - std::cout << "+++ Found candidates with 4, 5, 6 and 7 clusters:\t"; - for (int count : roadCounters) - std::cout << count << "\t"; - std::cout << std::endl; - - std::cout << "+++ Fitted candidates with 4, 5, 6 and 7 clusters:\t"; - for (int count : fitCounters) - std::cout << count << "\t"; - std::cout << std::endl; - - std::cout << "+++ Backprop candidates with 4, 5, 6 and 7 clusters:\t"; - for (int count : backpropagatedCounters) - std::cout << count << "\t"; - std::cout << std::endl; - - std::cout << "+++ Refitted candidates with 4, 5, 6 and 7 clusters:\t"; - for (int count : refitCounters) - std::cout << count << "\t"; - std::cout << std::endl; - - // std::cout << "+++ Cross check counters for 4, 5, 6 and 7 clusters:\t"; - // for (size_t iCount = 0; iCount < refitCounters.size(); ++iCount) { - // std::cout << xcheckCounters[iCount] << "\t"; - // //assert(refitCounters[iCount] == xcheckCounters[iCount]); - // } - // std::cout << std::endl; - - // std::cout << "+++ Nonsharing candidates with 4, 5, 6 and 7 clusters:\t"; - // for (int count : nonsharingCounters) - // std::cout << count << "\t"; - // std::cout << std::endl; - - // std::cout << "+++ Sharing matrix:\n"; - // for (int iCl = 4; iCl <= 7; ++iCl) { - // std::cout << "+++ "; - // for (int iSh = cumulativeIndex(iCl); iSh < cumulativeIndex(iCl + 1); ++iSh) { - // std::cout << sharingMatrix[iSh] << "\t"; - // } - // std::cout << std::endl; - // } -#endif } -bool Tracker::fitTrack(const ROframe& event, TrackITSExt& track, int start, int end, int step, const float chi2cut) +bool Tracker::fitTrack(TrackITSExt& track, int start, int end, int step, const float chi2cut, const float maxQoverPt) { auto propInstance = o2::base::Propagator::Instance(); track.setChi2(0); + int nCl{0}; for (int iLayer{start}; iLayer != end; iLayer += step) { if (track.getClusterIndex(iLayer) == constants::its::UnusedIndex) { continue; } - const TrackingFrameInfo& trackingHit = event.getTrackingFrameInfoOnLayer(iLayer).at(track.getClusterIndex(iLayer)); + const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer).at(track.getClusterIndex(iLayer)); if (!track.rotate(trackingHit.alphaTrackingFrame)) { return false; @@ -429,34 +370,35 @@ bool Tracker::fitTrack(const ROframe& event, TrackITSExt& track, int start, int } auto predChi2{track.getPredictedChi2(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)}; - if (predChi2 > chi2cut) { + if (nCl >= 3 && predChi2 > chi2cut * (nCl * 2 - 5)) { return false; } track.setChi2(track.getChi2() + predChi2); if (!track.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) { return false; } + nCl++; } - return true; + return std::abs(track.getQ2Pt()) < maxQoverPt; } void Tracker::traverseCellsTree(const int currentCellId, const int currentLayerId) { - Cell& currentCell{mPrimaryVertexContext->getCells()[currentLayerId][currentCellId]}; + Cell& currentCell{mTimeFrame->getCells()[currentLayerId][currentCellId]}; const int currentCellLevel = currentCell.getLevel(); - mPrimaryVertexContext->getRoads().back().addCell(currentLayerId, currentCellId); + mTimeFrame->getRoads().back().addCell(currentLayerId, currentCellId); if (currentLayerId > 0 && currentCellLevel > 1) { const int cellNeighboursNum{static_cast( - mPrimaryVertexContext->getCellsNeighbours()[currentLayerId - 1][currentCellId].size())}; + mTimeFrame->getCellsNeighbours()[currentLayerId - 1][currentCellId].size())}; bool isFirstValidNeighbour = true; for (int iNeighbourCell{0}; iNeighbourCell < cellNeighboursNum; ++iNeighbourCell) { const int neighbourCellId = - mPrimaryVertexContext->getCellsNeighbours()[currentLayerId - 1][currentCellId][iNeighbourCell]; - const Cell& neighbourCell = mPrimaryVertexContext->getCells()[currentLayerId - 1][neighbourCellId]; + mTimeFrame->getCellsNeighbours()[currentLayerId - 1][currentCellId][iNeighbourCell]; + const Cell& neighbourCell = mTimeFrame->getCells()[currentLayerId - 1][neighbourCellId]; if (currentCellLevel - 1 != neighbourCell.getLevel()) { continue; @@ -465,7 +407,7 @@ void Tracker::traverseCellsTree(const int currentCellId, const int currentLayerI if (isFirstValidNeighbour) { isFirstValidNeighbour = false; } else { - mPrimaryVertexContext->getRoads().push_back(mPrimaryVertexContext->getRoads().back()); + mTimeFrame->getRoads().push_back(mTimeFrame->getRoads().back()); } traverseCellsTree(neighbourCellId, currentLayerId - 1); @@ -476,20 +418,20 @@ void Tracker::traverseCellsTree(const int currentCellId, const int currentLayerI // currentCell.setLevel(0); } -void Tracker::computeRoadsMClabels(const ROframe& event) +void Tracker::computeRoadsMClabels() { /// Moore's Voting Algorithm - if (!event.hasMCinformation()) { + if (!mTimeFrame->hasMCinformation()) { return; } - mPrimaryVertexContext->initialiseRoadLabels(); + mTimeFrame->initialiseRoadLabels(); - int roadsNum{static_cast(mPrimaryVertexContext->getRoads().size())}; + int roadsNum{static_cast(mTimeFrame->getRoads().size())}; for (int iRoad{0}; iRoad < roadsNum; ++iRoad) { - Road& currentRoad{mPrimaryVertexContext->getRoads()[iRoad]}; + Road& currentRoad{mTimeFrame->getRoads()[iRoad]}; MCCompLabel maxOccurrencesValue{constants::its::UnusedIndex, constants::its::UnusedIndex, constants::its::UnusedIndex, false}; int count{0}; @@ -507,17 +449,18 @@ void Tracker::computeRoadsMClabels(const ROframe& event) } } - const Cell& currentCell{mPrimaryVertexContext->getCells()[iCell][currentCellIndex]}; + const Cell& currentCell{mTimeFrame->getCells()[iCell][currentCellIndex]}; if (isFirstRoadCell) { - const int cl0index{mPrimaryVertexContext->getClusters()[iCell][currentCell.getFirstClusterIndex()].clusterId}; - auto& cl0labs{event.getClusterLabels(iCell, cl0index)}; + const int cl0index{mTimeFrame->getClusters()[iCell][currentCell.getFirstClusterIndex()].clusterId}; + auto& cl0labs{mTimeFrame->getClusterLabels(iCell, cl0index)}; maxOccurrencesValue = cl0labs; count = 1; - const int cl1index{mPrimaryVertexContext->getClusters()[iCell + 1][currentCell.getSecondClusterIndex()].clusterId}; - const auto& cl1labs{event.getClusterLabels(iCell + 1, cl1index)}; + const int cl1index{mTimeFrame->getClusters()[iCell + 1][currentCell.getSecondClusterIndex()].clusterId}; + auto& cl1labs{mTimeFrame->getClusterLabels(iCell + 1, cl1index)}; + const int secondMonteCarlo{cl1labs.getTrackID()}; if (cl1labs == maxOccurrencesValue) { ++count; @@ -530,8 +473,9 @@ void Tracker::computeRoadsMClabels(const ROframe& event) isFirstRoadCell = false; } - const int cl2index{mPrimaryVertexContext->getClusters()[iCell + 2][currentCell.getThirdClusterIndex()].clusterId}; - const auto& cl2labs{event.getClusterLabels(iCell + 2, cl2index)}; + const int cl2index{mTimeFrame->getClusters()[iCell + 2][currentCell.getThirdClusterIndex()].clusterId}; + auto& cl2labs{mTimeFrame->getClusterLabels(iCell + 2, cl2index)}; + const int currentMonteCarlo = {cl2labs.getTrackID()}; if (cl2labs == maxOccurrencesValue) { ++count; @@ -546,62 +490,69 @@ void Tracker::computeRoadsMClabels(const ROframe& event) } } - mPrimaryVertexContext->setRoadLabel(iRoad, maxOccurrencesValue.getRawValue(), isFakeRoad); + mTimeFrame->setRoadLabel(iRoad, maxOccurrencesValue.getRawValue(), isFakeRoad); } } -void Tracker::computeTracksMClabels(const ROframe& event) +void Tracker::computeTracksMClabels() { /// Moore's Voting Algorithm - if (!event.hasMCinformation()) { + if (!mTimeFrame->hasMCinformation()) { return; } - int tracksNum{static_cast(mTracks.size())}; - - for (auto& track : mTracks) { - - MCCompLabel maxOccurrencesValue{constants::its::UnusedIndex, constants::its::UnusedIndex, - constants::its::UnusedIndex, false}; - int count{0}; - bool isFakeTrack{false}; + for (int iROF{0}; iROF < mTimeFrame->getNrof(); ++iROF) { + for (auto& track : mTimeFrame->getTracks(iROF)) { + MCCompLabel maxOccurrencesValue{constants::its::UnusedIndex, constants::its::UnusedIndex, + constants::its::UnusedIndex, false}; + int count{0}; + bool isFakeTrack{false}; + for (int iCluster = 0; iCluster < TrackITSExt::MaxClusters; ++iCluster) { + const int index = track.getClusterIndex(iCluster); + if (index == constants::its::UnusedIndex) { + continue; + } + const MCCompLabel& currentLabel = mTimeFrame->getClusterLabels(iCluster, index); + if (currentLabel == maxOccurrencesValue) { + ++count; + } else { + if (count != 0) { // only in the first iteration count can be 0 at this point + isFakeTrack = true; + --count; + } - for (int iCluster = 0; iCluster < TrackITSExt::MaxClusters; ++iCluster) { - const int index = track.getClusterIndex(iCluster); - if (index == constants::its::UnusedIndex) { - continue; - } - const MCCompLabel& currentLabel = event.getClusterLabels(iCluster, index); - if (currentLabel == maxOccurrencesValue) { - ++count; - } else { - if (count != 0) { // only in the first iteration count can be 0 at this point - isFakeTrack = true; - --count; + if (currentLabel == maxOccurrencesValue) { + ++count; + } else { + if (count != 0) { // only in the first iteration count can be 0 at this point + isFakeTrack = true; + --count; + } + if (count == 0) { + maxOccurrencesValue = currentLabel; + count = 1; + } + } } - if (count == 0) { - maxOccurrencesValue = currentLabel; - count = 1; + + if (isFakeTrack) { + maxOccurrencesValue.setFakeFlag(); } } - track.setExternalClusterIndex(iCluster, event.getClusterExternalIndex(iCluster, index)); - } - - if (isFakeTrack) { - maxOccurrencesValue.setFakeFlag(); + mTimeFrame->getTracksLabel(iROF).emplace_back(maxOccurrencesValue); } - mTrackLabels.emplace_back(maxOccurrencesValue); } } -void Tracker::rectifyClusterIndices(const ROframe& event) +void Tracker::rectifyClusterIndices() { - int tracksNum{static_cast(mTracks.size())}; - for (auto& track : mTracks) { - for (int iCluster = 0; iCluster < TrackITSExt::MaxClusters; ++iCluster) { - const int index = track.getClusterIndex(iCluster); - if (index != constants::its::UnusedIndex) { - track.setExternalClusterIndex(iCluster, event.getClusterExternalIndex(iCluster, index)); + for (int iROF{0}; iROF < mTimeFrame->getNrof(); ++iROF) { + for (auto& track : mTimeFrame->getTracks(iROF)) { + for (int iCluster = 0; iCluster < TrackITSExt::MaxClusters; ++iCluster) { + const int index = track.getClusterIndex(iCluster); + if (index != constants::its::UnusedIndex) { + track.setExternalClusterIndex(iCluster, mTimeFrame->getClusterExternalIndex(iCluster, index)); + } } } } @@ -630,7 +581,7 @@ track::TrackParCov Tracker::buildTrackSeed(const Cluster& cluster1, const Cluste const float tgl12 = math_utils::computeTanDipAngle(x1, y1, x2, y2, z1, z2); const float tgl23 = math_utils::computeTanDipAngle(x2, y2, x3, y3, z2, z3); - const float fy = 1. / (cluster2.rCoordinate - cluster3.rCoordinate); + const float fy = 1. / (cluster2.radius - cluster3.radius); const float& tz = fy; const float cy = (math_utils::computeCurvature(x1, y1, x2, y2 + constants::its::Resolution, x3, y3) - crv) / (constants::its::Resolution * getBz() * o2::constants::math::B2C) * diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraitsCPU.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraitsCPU.cxx index db1897dfd927b..cb8df35162e52 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraitsCPU.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraitsCPU.cxx @@ -32,240 +32,303 @@ namespace o2 namespace its { +constexpr int debugLevel{0}; + void TrackerTraitsCPU::computeLayerTracklets() { - PrimaryVertexContext* primaryVertexContext = mPrimaryVertexContext; - for (int iLayer{0}; iLayer < mTrkParams.TrackletsPerRoad(); ++iLayer) { - if (primaryVertexContext->getClusters()[iLayer].empty() || primaryVertexContext->getClusters()[iLayer + 1].empty()) { - continue; - } + TimeFrame* tf = mTimeFrame; - const float3& primaryVertex = primaryVertexContext->getPrimaryVertex(); - const int currentLayerClustersNum{static_cast(primaryVertexContext->getClusters()[iLayer].size())}; - - for (int iCluster{0}; iCluster < currentLayerClustersNum; ++iCluster) { - const Cluster& currentCluster{primaryVertexContext->getClusters()[iLayer][iCluster]}; - - if (primaryVertexContext->isClusterUsed(iLayer, currentCluster.clusterId)) { + for (int rof0{0}; rof0 < tf->getNrof(); ++rof0) { + gsl::span primaryVertices = tf->getPrimaryVertices(rof0); + for (int iLayer{0}; iLayer < mTrkParams.TrackletsPerRoad(); ++iLayer) { + gsl::span layer0 = tf->getClustersOnLayer(rof0, iLayer); + if (layer0.empty()) { continue; } - const float tanLambda{(currentCluster.zCoordinate - primaryVertex.z) / currentCluster.rCoordinate}; - const float zAtRmin{tanLambda * (mPrimaryVertexContext->getMinR(iLayer + 1) - - currentCluster.rCoordinate) + - currentCluster.zCoordinate}; - const float zAtRmax{tanLambda * (mPrimaryVertexContext->getMaxR(iLayer + 1) - - currentCluster.rCoordinate) + - currentCluster.zCoordinate}; - - const int4 selectedBinsRect{getBinsRect(currentCluster, iLayer, zAtRmin, zAtRmax, - mTrkParams.TrackletMaxDeltaZ[iLayer], mTrkParams.TrackletMaxDeltaPhi)}; - - if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) { - continue; - } - - int phiBinsNum{selectedBinsRect.w - selectedBinsRect.y + 1}; - - if (phiBinsNum < 0) { - phiBinsNum += mTrkParams.PhiBins; - } + const int currentLayerClustersNum{static_cast(layer0.size())}; + for (int iCluster{0}; iCluster < currentLayerClustersNum; ++iCluster) { + const Cluster& currentCluster{layer0[iCluster]}; + const int currentSortedIndex{tf->getSortedIndex(rof0, iLayer, iCluster)}; - for (int iPhiBin{selectedBinsRect.y}, iPhiCount{0}; iPhiCount < phiBinsNum; - iPhiBin = ++iPhiBin == mTrkParams.PhiBins ? 0 : iPhiBin, iPhiCount++) { - const int firstBinIndex{primaryVertexContext->mIndexTableUtils.getBinIndex(selectedBinsRect.x, iPhiBin)}; - const int maxBinIndex{firstBinIndex + selectedBinsRect.z - selectedBinsRect.x + 1}; - const int firstRowClusterIndex = primaryVertexContext->getIndexTables()[iLayer][firstBinIndex]; - const int maxRowClusterIndex = primaryVertexContext->getIndexTables()[iLayer][maxBinIndex]; + if (tf->isClusterUsed(iLayer, currentCluster.clusterId)) { + continue; + } - for (int iNextLayerCluster{firstRowClusterIndex}; iNextLayerCluster < maxRowClusterIndex; - ++iNextLayerCluster) { + for (auto& primaryVertex : primaryVertices) { + const float tanLambda{(currentCluster.zCoordinate - primaryVertex.getZ()) / currentCluster.radius}; - if (iNextLayerCluster >= (int)primaryVertexContext->getClusters()[iLayer + 1].size()) { - break; - } + const float zAtRmin{tanLambda * (tf->getMinR(iLayer + 1) - currentCluster.radius) + currentCluster.zCoordinate}; + const float zAtRmax{tanLambda * (tf->getMaxR(iLayer + 1) - currentCluster.radius) + currentCluster.zCoordinate}; - const Cluster& nextCluster{primaryVertexContext->getClusters()[iLayer + 1][iNextLayerCluster]}; + const int4 selectedBinsRect{getBinsRect(currentCluster, iLayer, zAtRmin, zAtRmax, + mTrkParams.TrackletMaxDeltaZ[iLayer], mTrkParams.TrackletMaxDeltaPhi)}; - if (primaryVertexContext->isClusterUsed(iLayer + 1, nextCluster.clusterId)) { + if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) { continue; } - const float deltaZ{o2::gpu::GPUCommonMath::Abs(tanLambda * (nextCluster.rCoordinate - currentCluster.rCoordinate) + - currentCluster.zCoordinate - nextCluster.zCoordinate)}; - const float deltaPhi{o2::gpu::GPUCommonMath::Abs(currentCluster.phiCoordinate - nextCluster.phiCoordinate)}; + int phiBinsNum{selectedBinsRect.w - selectedBinsRect.y + 1}; - if (deltaZ < mTrkParams.TrackletMaxDeltaZ[iLayer] && - (deltaPhi < mTrkParams.TrackletMaxDeltaPhi || - o2::gpu::GPUCommonMath::Abs(deltaPhi - constants::math::TwoPi) < mTrkParams.TrackletMaxDeltaPhi)) { - - if (iLayer > 0 && - primaryVertexContext->getTrackletsLookupTable()[iLayer - 1][iCluster] == constants::its::UnusedIndex) { + if (phiBinsNum < 0) { + phiBinsNum += mTrkParams.PhiBins; + } - primaryVertexContext->getTrackletsLookupTable()[iLayer - 1][iCluster] = - primaryVertexContext->getTracklets()[iLayer].size(); + int minRof = (rof0 >= mTrkParams.DeltaROF) ? rof0 - mTrkParams.DeltaROF : 0; + int maxRof = (rof0 == tf->getNrof() - mTrkParams.DeltaROF) ? rof0 : rof0 + mTrkParams.DeltaROF; + for (int rof1{minRof}; rof1 <= maxRof; ++rof1) { + gsl::span layer1 = tf->getClustersOnLayer(rof1, iLayer + 1); + if (layer1.empty()) { + continue; } - primaryVertexContext->getTracklets()[iLayer].emplace_back(iCluster, iNextLayerCluster, currentCluster, - nextCluster); + for (int iPhiBin{selectedBinsRect.y}, iPhiCount{0}; iPhiCount < phiBinsNum; + iPhiBin = (++iPhiBin == tf->mIndexTableUtils.getNphiBins()) ? 0 : iPhiBin, iPhiCount++) { + const int firstBinIndex{tf->mIndexTableUtils.getBinIndex(selectedBinsRect.x, iPhiBin)}; + const int maxBinIndex{firstBinIndex + selectedBinsRect.z - selectedBinsRect.x + 1}; + if constexpr (debugLevel) { + if (firstBinIndex < 0 || firstBinIndex > tf->getIndexTables(rof1)[iLayer].size() || + maxBinIndex < 0 || maxBinIndex > tf->getIndexTables(rof1)[iLayer].size()) { + std::cout << iLayer << "\t" << iCluster << "\t" << zAtRmin << "\t" << zAtRmax << "\t" << mTrkParams.TrackletMaxDeltaZ[iLayer] << "\t" << mTrkParams.TrackletMaxDeltaPhi << std::endl; + std::cout << currentCluster.zCoordinate << "\t" << primaryVertex.getZ() << "\t" << currentCluster.radius << std::endl; + std::cout << tf->getMinR(iLayer + 1) << "\t" << currentCluster.radius << "\t" << currentCluster.zCoordinate << std::endl; + std::cout << "Illegal access to IndexTable " << firstBinIndex << "\t" << maxBinIndex << "\t" << selectedBinsRect.z << "\t" << selectedBinsRect.x << std::endl; + exit(1); + } + } + const int firstRowClusterIndex = tf->getIndexTables(rof1)[iLayer][firstBinIndex]; + const int maxRowClusterIndex = tf->getIndexTables(rof1)[iLayer][maxBinIndex]; + + for (int iNextCluster{firstRowClusterIndex}; iNextCluster < maxRowClusterIndex; ++iNextCluster) { + if (iNextCluster >= (int)tf->getClusters()[iLayer + 1].size()) { + break; + } + const Cluster& nextCluster{layer1[iNextCluster]}; + + if (tf->isClusterUsed(iLayer + 1, nextCluster.clusterId)) { + continue; + } + + const float deltaZ{gpu::GPUCommonMath::Abs(tanLambda * (nextCluster.radius - currentCluster.radius) + + currentCluster.zCoordinate - nextCluster.zCoordinate)}; + const float deltaPhi{gpu::GPUCommonMath::Abs(currentCluster.phi - nextCluster.phi)}; + + if (deltaZ < mTrkParams.TrackletMaxDeltaZ[iLayer] && + (deltaPhi < mTrkParams.TrackletMaxDeltaPhi || + gpu::GPUCommonMath::Abs(deltaPhi - constants::math::TwoPi) < mTrkParams.TrackletMaxDeltaPhi)) { + if (iLayer > 0) { + tf->getTrackletsLookupTable()[iLayer - 1][currentSortedIndex]++; + } + tf->getTracklets()[iLayer].emplace_back(currentSortedIndex, tf->getSortedIndex(rof1, iLayer + 1, iNextCluster), currentCluster, + nextCluster, rof0, rof1); + } + } + } } } } } - if (iLayer > 0 && iLayer < mTrkParams.TrackletsPerRoad() - 1 && - primaryVertexContext->getTracklets()[iLayer].size() > primaryVertexContext->getCellsLookupTable()[iLayer - 1].size()) { - throw std::runtime_error(fmt::format("not enough memory in the CellsLookupTable, increase the tracklet memory coefficients: {} tracklets on L{}, lookup table size {} on L{}", - primaryVertexContext->getTracklets()[iLayer].size(), iLayer, primaryVertexContext->getCellsLookupTable()[iLayer - 1].size(), iLayer - 1)); + } + /// Cold code, fixups + + for (int iLayer{0}; iLayer < mTrkParams.CellsPerRoad(); ++iLayer) { + /// Sort tracklets + auto& trkl{tf->getTracklets()[iLayer + 1]}; + std::sort(trkl.begin(), trkl.end(), [](const Tracklet& a, const Tracklet& b) { + return a.firstClusterIndex < b.firstClusterIndex || (a.firstClusterIndex == b.firstClusterIndex && a.secondClusterIndex < b.secondClusterIndex); + }); + /// Remove duplicates + auto& lut{tf->getTrackletsLookupTable()[iLayer]}; + int id0{-1}, id1{-1}; + std::vector newTrk; + newTrk.reserve(trkl.size()); + for (auto& trk : trkl) { + if (trk.firstClusterIndex == id0 && trk.secondClusterIndex == id1) { + lut[id0]--; + } else { + id0 = trk.firstClusterIndex; + id1 = trk.secondClusterIndex; + newTrk.push_back(trk); + } + } + trkl.swap(newTrk); + + /// Compute LUT + std::exclusive_scan(lut.begin(), lut.end(), lut.begin(), 0); + lut.push_back(trkl.size()); + } + /// Layer 0 is done outside the loop + std::sort(tf->getTracklets()[0].begin(), tf->getTracklets()[0].end(), [](const Tracklet& a, const Tracklet& b) { + return a.firstClusterIndex < b.firstClusterIndex || (a.firstClusterIndex == b.firstClusterIndex && a.secondClusterIndex < b.secondClusterIndex); + }); + int id0{-1}, id1{-1}; + std::vector newTrk; + newTrk.reserve(tf->getTracklets()[0].size()); + for (auto& trk : tf->getTracklets()[0]) { + if (trk.firstClusterIndex != id0 || trk.secondClusterIndex != id1) { + id0 = trk.firstClusterIndex; + id1 = trk.secondClusterIndex; + newTrk.push_back(trk); } } -#ifdef CA_DEBUG - std::cout << "+++ Number of tracklets per layer: "; - for (int iLayer{0}; iLayer < mTrkParams.TrackletsPerRoad(); ++iLayer) { - std::cout << primaryVertexContext->getTracklets()[iLayer].size() << "\t"; + tf->getTracklets()[0].swap(newTrk); + + /// Create tracklets labels + if (tf->hasMCinformation()) { + for (int iLayer{0}; iLayer < mTrkParams.TrackletsPerRoad(); ++iLayer) { + for (auto& trk : tf->getTracklets()[iLayer]) { + int currentId{tf->getClusters()[iLayer][trk.firstClusterIndex].clusterId}; + int nextId{tf->getClusters()[iLayer + 1][trk.secondClusterIndex].clusterId}; + MCCompLabel currentLab{tf->getClusterLabels(iLayer, currentId)}; + MCCompLabel nextLab{tf->getClusterLabels(iLayer + 1, nextId)}; + tf->getTrackletsLabel(iLayer).emplace_back(currentLab == nextLab ? currentLab : MCCompLabel()); + } + } } - std::cout << std::endl; -#endif } void TrackerTraitsCPU::computeLayerCells() { - PrimaryVertexContext* primaryVertexContext = mPrimaryVertexContext; + TimeFrame* tf = mTimeFrame; for (int iLayer{0}; iLayer < mTrkParams.CellsPerRoad(); ++iLayer) { - if (primaryVertexContext->getTracklets()[iLayer + 1].empty() || - primaryVertexContext->getTracklets()[iLayer].empty()) { - - return; + if (tf->getTracklets()[iLayer + 1].empty() || + tf->getTracklets()[iLayer].empty()) { + continue; } - const float3& primaryVertex = primaryVertexContext->getPrimaryVertex(); - const int currentLayerTrackletsNum{static_cast(primaryVertexContext->getTracklets()[iLayer].size())}; + const int currentLayerTrackletsNum{static_cast(tf->getTracklets()[iLayer].size())}; for (int iTracklet{0}; iTracklet < currentLayerTrackletsNum; ++iTracklet) { - const Tracklet& currentTracklet{primaryVertexContext->getTracklets()[iLayer][iTracklet]}; + const Tracklet& currentTracklet{tf->getTracklets()[iLayer][iTracklet]}; const int nextLayerClusterIndex{currentTracklet.secondClusterIndex}; const int nextLayerFirstTrackletIndex{ - primaryVertexContext->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex]}; - - if (nextLayerFirstTrackletIndex == constants::its::UnusedIndex) { + tf->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex]}; + const int nextLayerLastTrackletIndex{ + tf->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex + 1]}; + if (nextLayerFirstTrackletIndex == nextLayerLastTrackletIndex) { continue; } - const Cluster& firstCellCluster{primaryVertexContext->getClusters()[iLayer][currentTracklet.firstClusterIndex]}; - const Cluster& secondCellCluster{ - primaryVertexContext->getClusters()[iLayer + 1][currentTracklet.secondClusterIndex]}; - const float firstCellClusterQuadraticRCoordinate{firstCellCluster.rCoordinate * firstCellCluster.rCoordinate}; - const float secondCellClusterQuadraticRCoordinate{secondCellCluster.rCoordinate * - secondCellCluster.rCoordinate}; - const float3 firstDeltaVector{secondCellCluster.xCoordinate - firstCellCluster.xCoordinate, - secondCellCluster.yCoordinate - firstCellCluster.yCoordinate, - secondCellClusterQuadraticRCoordinate - firstCellClusterQuadraticRCoordinate}; - const int nextLayerTrackletsNum{static_cast(primaryVertexContext->getTracklets()[iLayer + 1].size())}; - - for (int iNextLayerTracklet{nextLayerFirstTrackletIndex}; - iNextLayerTracklet < nextLayerTrackletsNum && - primaryVertexContext->getTracklets()[iLayer + 1][iNextLayerTracklet].firstClusterIndex == - nextLayerClusterIndex; - ++iNextLayerTracklet) { - - const Tracklet& nextTracklet{primaryVertexContext->getTracklets()[iLayer + 1][iNextLayerTracklet]}; + const Cluster& cellClus0{tf->getClusters()[iLayer][currentTracklet.firstClusterIndex]}; + const Cluster& cellClus1{ + tf->getClusters()[iLayer + 1][currentTracklet.secondClusterIndex]}; + const float cellClus0R2{cellClus0.radius * cellClus0.radius}; + const float cellClus1R2{cellClus1.radius * cellClus1.radius}; + const float3 firstDeltaVector{cellClus1.xCoordinate - cellClus0.xCoordinate, + cellClus1.yCoordinate - cellClus0.yCoordinate, + cellClus1R2 - cellClus0R2}; + + for (int iNextTracklet{nextLayerFirstTrackletIndex}; iNextTracklet < nextLayerLastTrackletIndex; ++iNextTracklet) { + if (tf->getTracklets()[iLayer + 1][iNextTracklet].firstClusterIndex != nextLayerClusterIndex) { + break; + } + const Tracklet& nextTracklet{tf->getTracklets()[iLayer + 1][iNextTracklet]}; const float deltaTanLambda{std::abs(currentTracklet.tanLambda - nextTracklet.tanLambda)}; - const float deltaPhi{std::abs(currentTracklet.phiCoordinate - nextTracklet.phiCoordinate)}; + const float deltaPhi{std::abs(currentTracklet.phi - nextTracklet.phi)}; if (deltaTanLambda < mTrkParams.CellMaxDeltaTanLambda && (deltaPhi < mTrkParams.CellMaxDeltaPhi || std::abs(deltaPhi - constants::math::TwoPi) < mTrkParams.CellMaxDeltaPhi)) { const float averageTanLambda{0.5f * (currentTracklet.tanLambda + nextTracklet.tanLambda)}; - const float directionZIntersection{-averageTanLambda * firstCellCluster.rCoordinate + - firstCellCluster.zCoordinate}; - const float deltaZ{std::abs(directionZIntersection - primaryVertex.z)}; + const float directionZIntersection{-averageTanLambda * cellClus0.radius + + cellClus0.zCoordinate}; + + unsigned short romin = std::min(std::min(currentTracklet.rof[0], currentTracklet.rof[1]), nextTracklet.rof[1]); + unsigned short romax = std::max(std::max(currentTracklet.rof[0], currentTracklet.rof[1]), nextTracklet.rof[1]); + bool deltaZflag{false}; + gsl::span primaryVertices{tf->getPrimaryVertices(romin, romax)}; + for (const auto& primaryVertex : primaryVertices) { + deltaZflag |= std::abs(directionZIntersection - primaryVertex.getZ()) < mTrkParams.CellMaxDeltaZ[iLayer]; + } - if (deltaZ < mTrkParams.CellMaxDeltaZ[iLayer]) { + if (deltaZflag) { const Cluster& thirdCellCluster{ - primaryVertexContext->getClusters()[iLayer + 2][nextTracklet.secondClusterIndex]}; + tf->getClusters()[iLayer + 2][nextTracklet.secondClusterIndex]}; - const float thirdCellClusterQuadraticRCoordinate{thirdCellCluster.rCoordinate * - thirdCellCluster.rCoordinate}; + const float thirdCellClusterR2{thirdCellCluster.radius * + thirdCellCluster.radius}; - const float3 secondDeltaVector{thirdCellCluster.xCoordinate - firstCellCluster.xCoordinate, - thirdCellCluster.yCoordinate - firstCellCluster.yCoordinate, - thirdCellClusterQuadraticRCoordinate - - firstCellClusterQuadraticRCoordinate}; + const float3 secondDeltaVector{thirdCellCluster.xCoordinate - cellClus0.xCoordinate, + thirdCellCluster.yCoordinate - cellClus0.yCoordinate, + thirdCellClusterR2 - cellClus0R2}; float3 cellPlaneNormalVector{math_utils::crossProduct(firstDeltaVector, secondDeltaVector)}; - const float vectorNorm{std::sqrt(cellPlaneNormalVector.x * cellPlaneNormalVector.x + - cellPlaneNormalVector.y * cellPlaneNormalVector.y + - cellPlaneNormalVector.z * cellPlaneNormalVector.z)}; + const float vectorNorm{std::hypot(cellPlaneNormalVector.x, cellPlaneNormalVector.y, cellPlaneNormalVector.z)}; if (vectorNorm < constants::math::FloatMinThreshold || std::abs(cellPlaneNormalVector.z) < constants::math::FloatMinThreshold) { - continue; } const float inverseVectorNorm{1.0f / vectorNorm}; - const float3 normalizedPlaneVector{cellPlaneNormalVector.x * inverseVectorNorm, - cellPlaneNormalVector.y * inverseVectorNorm, - cellPlaneNormalVector.z * inverseVectorNorm}; - const float planeDistance{-normalizedPlaneVector.x * (secondCellCluster.xCoordinate - primaryVertex.x) - - (normalizedPlaneVector.y * secondCellCluster.yCoordinate - primaryVertex.y) - - normalizedPlaneVector.z * secondCellClusterQuadraticRCoordinate}; - const float normalizedPlaneVectorQuadraticZCoordinate{normalizedPlaneVector.z * normalizedPlaneVector.z}; - const float cellTrajectoryRadius{std::sqrt( - (1.0f - normalizedPlaneVectorQuadraticZCoordinate - 4.0f * planeDistance * normalizedPlaneVector.z) / - (4.0f * normalizedPlaneVectorQuadraticZCoordinate))}; - const float2 circleCenter{-0.5f * normalizedPlaneVector.x / normalizedPlaneVector.z, - -0.5f * normalizedPlaneVector.y / normalizedPlaneVector.z}; - const float distanceOfClosestApproach{std::abs( - cellTrajectoryRadius - std::sqrt(circleCenter.x * circleCenter.x + circleCenter.y * circleCenter.y))}; - - if (distanceOfClosestApproach > - mTrkParams.CellMaxDCA[iLayer]) { - + const float3 normVect{cellPlaneNormalVector.x * inverseVectorNorm, + cellPlaneNormalVector.y * inverseVectorNorm, + cellPlaneNormalVector.z * inverseVectorNorm}; + const float planeDistance{-normVect.x * (cellClus1.xCoordinate - tf->getBeamX()) - normVect.y * (cellClus1.yCoordinate - tf->getBeamY()) - normVect.z * cellClus1R2}; + const float normVectZsquare{normVect.z * normVect.z}; + const float cellRadius{std::sqrt( + (1.0f - normVectZsquare - 4.0f * planeDistance * normVect.z) / + (4.0f * normVectZsquare))}; + const float2 circleCenter{-0.5f * normVect.x / normVect.z, + -0.5f * normVect.y / normVect.z}; + const float dca{std::abs(cellRadius - std::hypot(circleCenter.x, circleCenter.y))}; + + if (dca > mTrkParams.CellMaxDCA[iLayer]) { continue; } - const float cellTrajectoryCurvature{1.0f / cellTrajectoryRadius}; - if (iLayer > 0 && - primaryVertexContext->getCellsLookupTable()[iLayer - 1][iTracklet] == constants::its::UnusedIndex) { - - primaryVertexContext->getCellsLookupTable()[iLayer - 1][iTracklet] = - primaryVertexContext->getCells()[iLayer].size(); + const float cellTrajectoryCurvature{1.0f / cellRadius}; + if (iLayer > 0 && tf->getCellsLookupTable()[iLayer - 1].size() <= iTracklet) { + tf->getCellsLookupTable()[iLayer - 1].resize(iTracklet + 1, tf->getCells()[iLayer].size()); } - primaryVertexContext->getCells()[iLayer].emplace_back( + tf->getCells()[iLayer].emplace_back( currentTracklet.firstClusterIndex, nextTracklet.firstClusterIndex, nextTracklet.secondClusterIndex, - iTracklet, iNextLayerTracklet, normalizedPlaneVector, cellTrajectoryCurvature); + iTracklet, iNextTracklet, normVect, cellTrajectoryCurvature); } } } } + if (iLayer > 0) { + tf->getCellsLookupTable()[iLayer - 1].resize(currentLayerTrackletsNum + 1, currentLayerTrackletsNum); + } } -#ifdef CA_DEBUG - std::cout << "+++ Number of cells per layer: "; - for (int iLayer{0}; iLayer < mTrkParams.CellsPerRoad(); ++iLayer) { - std::cout << primaryVertexContext->getCells()[iLayer].size() << "\t"; + + /// Create cells labels + if (tf->hasMCinformation()) { + for (int iLayer{0}; iLayer < mTrkParams.CellsPerRoad(); ++iLayer) { + for (auto& cell : tf->getCells()[iLayer]) { + MCCompLabel currentLab{tf->getTrackletsLabel(iLayer)[cell.getFirstTrackletIndex()]}; + MCCompLabel nextLab{tf->getTrackletsLabel(iLayer + 1)[cell.getSecondTrackletIndex()]}; + tf->getCellsLabel(iLayer).emplace_back(currentLab == nextLab ? currentLab : MCCompLabel()); + } + } + } + + if constexpr (debugLevel) { + for (int iLayer{0}; iLayer < mTrkParams.CellsPerRoad(); ++iLayer) { + std::cout << "Cells on layer " << iLayer << " " << tf->getCells()[iLayer].size() << std::endl; + } } - std::cout << std::endl; -#endif } void TrackerTraitsCPU::refitTracks(const std::vector>& tf, std::vector& tracks) { std::vector cells; for (int iLayer = 0; iLayer < mTrkParams.CellsPerRoad(); iLayer++) { - cells.push_back(mPrimaryVertexContext->getCells()[iLayer].data()); + cells.push_back(mTimeFrame->getCells()[iLayer].data()); } std::vector clusters; for (int iLayer = 0; iLayer < mTrkParams.NLayers; iLayer++) { - clusters.push_back(mPrimaryVertexContext->getClusters()[iLayer].data()); + clusters.push_back(mTimeFrame->getClusters()[iLayer].data()); } - mChainRunITSTrackFit(*mChain, mPrimaryVertexContext->getRoads(), clusters, cells, tf, tracks); + mChainRunITSTrackFit(*mChain, mTimeFrame->getRoads(), clusters, cells, tf, tracks); } } // namespace its diff --git a/Detectors/ITSMFT/ITS/tracking/src/Vertexer.cxx b/Detectors/ITSMFT/ITS/tracking/src/Vertexer.cxx index f8c723aa2f9ec..6c16ce4f40e13 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Vertexer.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Vertexer.cxx @@ -33,19 +33,19 @@ Vertexer::Vertexer(VertexerTraits* traits) mTraits = traits; } -float Vertexer::clustersToVertices(ROframe& event, const bool useMc, std::ostream& timeBenchmarkOutputStream) +float Vertexer::clustersToVertices(ROframe& event, const bool useMc, std::function logger) { ROframe* eventptr = &event; float total{0.f}; - total += evaluateTask(&Vertexer::initialiseVertexer, "Vertexer initialisation", timeBenchmarkOutputStream, eventptr); - total += evaluateTask(&Vertexer::findTracklets, "Tracklet finding", timeBenchmarkOutputStream); + total += evaluateTask(&Vertexer::initialiseVertexer, "Vertexer initialisation", logger, eventptr); + total += evaluateTask(&Vertexer::findTracklets, "Tracklet finding", logger); #ifdef _ALLOW_DEBUG_TREES_ITS_ if (useMc) { - total += evaluateTask(&Vertexer::filterMCTracklets, "MC tracklets filtering", timeBenchmarkOutputStream); + total += evaluateTask(&Vertexer::filterMCTracklets, "MC tracklets filtering", logger); } #endif - total += evaluateTask(&Vertexer::validateTracklets, "Adjacent tracklets validation", timeBenchmarkOutputStream); - total += evaluateTask(&Vertexer::findVertices, "Vertex finding", timeBenchmarkOutputStream); + total += evaluateTask(&Vertexer::validateTracklets, "Adjacent tracklets validation", logger); + total += evaluateTask(&Vertexer::findVertices, "Vertex finding", logger); return total; } diff --git a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx index 7b89545a07692..37db06873d17e 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx @@ -73,7 +73,7 @@ void trackleterKernelSerial( // loop on clusters next layer for (int iNextLayerClusterIndex{firstRowClusterIndex}; iNextLayerClusterIndex < maxRowClusterIndex && iNextLayerClusterIndex < static_cast(clustersNextLayer.size()); ++iNextLayerClusterIndex) { const Cluster& nextCluster{clustersNextLayer[iNextLayerClusterIndex]}; - if (o2::gpu::GPUCommonMath::Abs(currentCluster.phiCoordinate - nextCluster.phiCoordinate) < phiCut) { + if (o2::gpu::GPUCommonMath::Abs(currentCluster.phi - nextCluster.phi) < phiCut) { if (storedTracklets < maxTrackletsPerCluster) { if (pairOfLayers == LAYER0_TO_LAYER1) { Tracklets.emplace_back(iNextLayerClusterIndex, iCurrentLayerClusterIndex, nextCluster, currentCluster); @@ -114,7 +114,7 @@ void trackletSelectionKernelSerial( for (int iTracklet12{offset12}; iTracklet12 < offset12 + foundTracklets12[iCurrentLayerClusterIndex]; ++iTracklet12) { for (int iTracklet01{offset01}; iTracklet01 < offset01 + foundTracklets01[iCurrentLayerClusterIndex]; ++iTracklet01) { const float deltaTanLambda{o2::gpu::GPUCommonMath::Abs(tracklets01[iTracklet01].tanLambda - tracklets12[iTracklet12].tanLambda)}; - const float deltaPhi{o2::gpu::GPUCommonMath::Abs(tracklets01[iTracklet01].phiCoordinate - tracklets12[iTracklet12].phiCoordinate)}; + const float deltaPhi{o2::gpu::GPUCommonMath::Abs(tracklets01[iTracklet01].phi - tracklets12[iTracklet12].phi)}; if (deltaTanLambda < tanLambdaCut && deltaPhi < phiCut && validTracklets != maxTracklets) { assert(tracklets01[iTracklet01].secondClusterIndex == tracklets12[iTracklet12].firstClusterIndex); #ifdef _ALLOW_DEBUG_TREES_ITS_ @@ -199,7 +199,7 @@ void VertexerTraits::arrangeClusters(ROframe* event) } for (unsigned int iCluster{0}; iCluster < clustersNum; ++iCluster) { mClusters[iLayer].emplace_back(iLayer, mIndexTableUtils, currentLayer.at(iCluster)); - mAverageClustersRadii[iLayer] += mClusters[iLayer].back().rCoordinate; + mAverageClustersRadii[iLayer] += mClusters[iLayer].back().radius; } mAverageClustersRadii[iLayer] *= 1.f / clustersNum; diff --git a/Detectors/ITSMFT/ITS/workflow/include/ITSWorkflow/TrackerSpec.h b/Detectors/ITSMFT/ITS/workflow/include/ITSWorkflow/TrackerSpec.h index a43ac8edaafe2..a57e56df4dd47 100644 --- a/Detectors/ITSMFT/ITS/workflow/include/ITSWorkflow/TrackerSpec.h +++ b/Detectors/ITSMFT/ITS/workflow/include/ITSWorkflow/TrackerSpec.h @@ -20,6 +20,7 @@ #include "Framework/DataProcessorSpec.h" #include "Framework/Task.h" +#include "ITStracking/TimeFrame.h" #include "ITStracking/Tracker.h" #include "ITStracking/TrackerTraitsCPU.h" #include "ITStracking/Vertexer.h" @@ -50,6 +51,7 @@ class TrackerDPL : public framework::Task bool mRunVertexer = true; std::string mMode = "sync"; o2::itsmft::TopologyDictionary mDict; + TimeFrame mTimeFrame; std::unique_ptr mRecChain = nullptr; std::unique_ptr mGRP = nullptr; std::unique_ptr mTracker = nullptr; diff --git a/Detectors/ITSMFT/ITS/workflow/src/CookedTrackerSpec.cxx b/Detectors/ITSMFT/ITS/workflow/src/CookedTrackerSpec.cxx index b3e00b86a1451..6183d397347d4 100644 --- a/Detectors/ITSMFT/ITS/workflow/src/CookedTrackerSpec.cxx +++ b/Detectors/ITSMFT/ITS/workflow/src/CookedTrackerSpec.cxx @@ -153,7 +153,7 @@ void CookedTrackerDPL::run(ProcessingContext& pc) } } - vertexer.clustersToVertices(event); + vertexer.clustersToVertices(event, false, [&](std::string s) { LOG(INFO) << s; }); auto vtxVecLoc = vertexer.exportVertices(); if (multEstConf.cutMultVtxLow > 0 || multEstConf.cutMultVtxHigh > 0) { // cut was requested diff --git a/Detectors/ITSMFT/ITS/workflow/src/TrackerSpec.cxx b/Detectors/ITSMFT/ITS/workflow/src/TrackerSpec.cxx index daf00c5f91413..3379d254bbfd0 100644 --- a/Detectors/ITSMFT/ITS/workflow/src/TrackerSpec.cxx +++ b/Detectors/ITSMFT/ITS/workflow/src/TrackerSpec.cxx @@ -23,7 +23,7 @@ #include "SimulationDataFormat/MCCompLabel.h" #include "SimulationDataFormat/MCTruthContainer.h" #include "DataFormatsITSMFT/ROFRecord.h" -#include "CommonDataFormat/IRFrame.h" + #include "ITStracking/ROframe.h" #include "ITStracking/IOUtils.h" #include "ITStracking/TrackingConfigParam.h" @@ -33,6 +33,7 @@ #include "DetectorsBase/GeometryManager.h" #include "DetectorsBase/Propagator.h" #include "ITSBase/GeometryTGeo.h" +#include "CommonDataFormat/IRFrame.h" #include "DetectorsCommonDataFormats/NameConf.h" #include "ITSReconstruction/FastMultEstConfig.h" @@ -80,7 +81,7 @@ void TrackerDPL::init(InitContext& ic) auto* chainITS = mRecChain->AddChain(); mRecChain->Init(); mVertexer = std::make_unique(chainITS->GetITSVertexerTraits()); - mTracker = std::make_unique(chainITS->GetITSTrackerTraits()); + mTracker = std::make_unique(new TrackerTraitsCPU(&mTimeFrame)); std::vector trackParams; std::vector memParams; @@ -89,11 +90,14 @@ void TrackerDPL::init(InitContext& ic) if (mMode == "async") { trackParams.resize(3); - memParams.resize(3); trackParams[0].TrackletMaxDeltaPhi = 0.05f; - trackParams[1].TrackletMaxDeltaPhi = 0.1f; + 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; - trackParams[2].TrackletMaxDeltaPhi = 0.3; + memParams.resize(3); LOG(INFO) << "Initializing tracker in async. phase reconstruction with " << trackParams.size() << " passes"; } else if (mMode == "sync") { @@ -159,6 +163,11 @@ void TrackerDPL::run(ProcessingContext& pc) auto rofsinput = pc.inputs().get>("ROframes"); auto& rofs = pc.outputs().make>(Output{"ITS", "ITSTrackROF", 0, Lifetime::Timeframe}, rofsinput.begin(), rofsinput.end()); + auto& irFrames = pc.outputs().make>(Output{"ITS", "IRFRAMES", 0, Lifetime::Timeframe}); + + const auto& alpParams = o2::itsmft::DPLAlpideParam::Instance(); // RS: this should come from CCDB + int nBCPerTF = alpParams.roFrameLengthInBC; + LOG(INFO) << "ITSTracker pulled " << compClusters.size() << " clusters, " << rofs.size() << " RO frames"; const dataformats::MCTruthContainer* labels = nullptr; @@ -179,47 +188,28 @@ void TrackerDPL::run(ProcessingContext& pc) auto& vertROFvec = pc.outputs().make>(Output{"ITS", "VERTICESROF", 0, Lifetime::Timeframe}); auto& vertices = pc.outputs().make>(Output{"ITS", "VERTICES", 0, Lifetime::Timeframe}); - auto& irFrames = pc.outputs().make>(Output{"ITS", "IRFRAMES", 0, Lifetime::Timeframe}); - std::uint32_t roFrame = 0; ROframe event(0, 7); bool continuous = mGRP->isDetContinuousReadOut("ITS"); LOG(INFO) << "ITSTracker RO: continuous=" << continuous; - const auto& alpParams = o2::itsmft::DPLAlpideParam::Instance(); // RS: this should come from CCDB - int nBCPerTF = continuous ? alpParams.roFrameLengthInBC : alpParams.roFrameLengthTrig; - const auto& multEstConf = FastMultEstConfig::Instance(); // parameters for mult estimation and cuts FastMultEst multEst; // mult estimator - // snippet to convert found tracks to final output tracks with separate cluster indices - auto copyTracks = [](auto& tracks, auto& allTracks, auto& allClusIdx) { - for (auto& trc : tracks) { - 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); - } - }; - gsl::span::iterator pattIt = patterns.begin(); - for (auto& rof : rofs) { - int nclUsed = ioutils::loadROFrameData(rof, event, compClusters, pattIt, mDict, labels); - // prepare in advance output ROFRecords, even if this ROF to be rejected - int first = allTracks.size(); + if (continuous) { + gsl::span rofspan(rofs); + mTimeFrame.loadROFrameData(rofspan, compClusters, pattIt, mDict, labels); + pattIt = patterns.begin(); + std::vector savedROF; + auto logger = [&](std::string s) { LOG(INFO) << s; }; + for (auto& rof : rofspan) { + + int nclUsed = ioutils::loadROFrameData(rof, event, compClusters, pattIt, mDict, labels); + // prepare in advance output ROFRecords, even if this ROF to be rejected + int first = allTracks.size(); - if (nclUsed) { LOG(INFO) << "ROframe: " << roFrame << ", clusters loaded : " << nclUsed; // for vertices output @@ -227,70 +217,62 @@ void TrackerDPL::run(ProcessingContext& pc) vtxROF.setFirstEntry(vertices.size()); // dedicated ROFRecord vtxROF.setNEntries(0); - if (multEstConf.cutMultClusLow > 0 || multEstConf.cutMultClusHigh > 0) { // cut was requested - auto mult = multEst.process(rof.getROFData(compClusters)); - if (mult < multEstConf.cutMultClusLow || mult > multEstConf.cutMultClusHigh) { - LOG(INFO) << "Estimated cluster mult. " << mult << " is outside of requested range " - << multEstConf.cutMultClusLow << " : " << multEstConf.cutMultClusHigh << " | ROF " << rof.getBCData(); - rof.setFirstEntry(first); - rof.setNEntries(0); - continue; - } - } - std::vector vtxVecLoc; if (mRunVertexer) { - mVertexer->clustersToVertices(event); + mVertexer->clustersToVertices(event, false, logger); vtxVecLoc = mVertexer->exportVertices(); } + mTimeFrame.addPrimaryVertices(vtxVecLoc); - if (mRunVertexer && (multEstConf.cutMultVtxLow > 0 || multEstConf.cutMultVtxHigh > 0)) { // cut was requested - std::vector>> vtxVecSel; - vtxVecSel.swap(vtxVecLoc); - for (const auto& vtx : vtxVecSel) { - if (vtx.getNContributors() < multEstConf.cutMultVtxLow || (multEstConf.cutMultVtxHigh > 0 && vtx.getNContributors() > multEstConf.cutMultVtxHigh)) { - LOG(INFO) << "Found vertex mult. " << vtx.getNContributors() << " is outside of requested range " - << multEstConf.cutMultVtxLow << " : " << multEstConf.cutMultVtxHigh << " | ROF " << rof.getBCData(); - continue; // skip vertex of unwanted multiplicity - } - vtxVecLoc.push_back(vtx); - } - if (vtxVecLoc.empty()) { // reject ROF - rof.setFirstEntry(first); - rof.setNEntries(0); - continue; - } - } - - if (mRunVertexer) { - event.addPrimaryVertices(vtxVecLoc); - } else { - event.addPrimaryVertex(0.f, 0.f, 0.f); - } - mTracker->setROFrame(roFrame); - mTracker->clustersToTracks(event); - tracks.swap(mTracker->getTracks()); - LOG(INFO) << "Found tracks: " << tracks.size(); - int number = tracks.size(); - trackLabels.swap(mTracker->getTrackLabels()); /// FIXME: assignment ctor is not optimal. - rof.setFirstEntry(first); - rof.setNEntries(number); - copyTracks(tracks, allTracks, allClusIdx); - std::copy(trackLabels.begin(), trackLabels.end(), std::back_inserter(allTrackLabels)); - trackLabels.clear(); vtxROF.setNEntries(vtxVecLoc.size()); for (const auto& vtx : vtxVecLoc) { vertices.push_back(vtx); } - if (number) { + savedROF.push_back(roFrame); + roFrame++; + } + mTracker->clustersToTracks(logger); + + for (unsigned int iROF{0}; iROF < rofs.size(); ++iROF) { + + auto& rof{rofs[iROF]}; + tracks = mTimeFrame.getTracks(iROF); + trackLabels = mTimeFrame.getTracksLabel(iROF); + auto number{tracks.size()}; + auto first{allTracks.size()}; + int offset = -rof.getFirstEntry(); // cluster entry!!! + rof.setFirstEntry(first); + rof.setNEntries(number); + + if (tracks.size()) { irFrames.emplace_back(rof.getBCData(), rof.getBCData() + nBCPerTF - 1); } + std::copy(trackLabels.begin(), trackLabels.end(), std::back_inserter(allTrackLabels)); + // Some conversions that needs to be moved in the tracker internals + for (unsigned int iTrk{0}; iTrk < tracks.size(); ++iTrk) { + 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); + } } - roFrame++; - } + } //MP: no triggered mode for the time being LOG(INFO) << "ITSTracker pushed " << allTracks.size() << " tracks"; if (mIsMC) { + LOG(INFO) << "ITSTracker pushed " << allTrackLabels.size() << " track labels"; + pc.outputs().snapshot(Output{"ITS", "TRACKSMCTR", 0, Lifetime::Timeframe}, allTrackLabels); pc.outputs().snapshot(Output{"ITS", "ITSTrackMC2ROF", 0, Lifetime::Timeframe}, mc2rofs); } diff --git a/Detectors/Upgrades/ALICE3/TRK/macros/ALICE3toAO2D.C b/Detectors/Upgrades/ALICE3/TRK/macros/ALICE3toAO2D.C index f7ebbe8dc9129..971682b5f65b4 100644 --- a/Detectors/Upgrades/ALICE3/TRK/macros/ALICE3toAO2D.C +++ b/Detectors/Upgrades/ALICE3/TRK/macros/ALICE3toAO2D.C @@ -34,6 +34,7 @@ #include "ITStracking/IOUtils.h" #include "ITStracking/Tracker.h" #include "ITStracking/TrackerTraitsCPU.h" +#include "ITStracking/TimeFrame.h" #include "ITStracking/Vertexer.h" #include "ITStracking/VertexerTraits.h" #include "DataFormatsITSMFT/Cluster.h" @@ -99,7 +100,8 @@ void ALICE3toAO2D() itsHits.AddFile(hitsFileName.data()); - o2::its::Tracker tracker(new o2::its::TrackerTraitsCPU()); + o2::its::TimeFrame tf; + o2::its::Tracker tracker(new o2::its::TrackerTraitsCPU(&tf)); tracker.setBz(5.f); std::uint32_t roFrame; @@ -439,10 +441,16 @@ void ALICE3toAO2D() std::vector vertices = vertexer.exportVertices(); std::cout << "*- Number of vertices found: " << vertices.size() << endl; hNVertices->Fill(vertices.size()); - //Skip events with no vertex + + tf.addPrimaryVertices(vertices); + } + + tracker.clustersToTracks(); + + for (int iROF{0}; iROF < tf.getNrof(); ++iROF) { + auto vertices{tf.getPrimaryVertices(iROF)}; if (vertices.size() == 0) { std::cout << "*- No primary vertex found, skipping event" << std::endl; - continue; } o2::math_utils::Point3D pos{vertices[0].getX(), vertices[0].getY(), vertices[0].getZ()}; @@ -501,11 +509,8 @@ void ALICE3toAO2D() bc.fTriggerMask |= 1ull << iii; bc.fRunNumber = 246087; //ah, the good old days - std::cout << "*- Event " << iEvent << " tracking" << std::endl; - tracker.clustersToTracks(event); - auto& lTracks = tracker.getTracks(); - auto& lTracksLabels = tracker.getTrackLabels(); - std::cout << "*- Event " << iEvent << " done tracking!" << std::endl; + auto& lTracks = tf.getTracks(iROF); + auto& lTracksLabels = tf.getTracksLabel(iROF); hNTracks->Fill(lTracks.size()); //+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ diff --git a/Detectors/Upgrades/IT3/README.md b/Detectors/Upgrades/IT3/README.md index 9a4d82a533edb..f8643db684fb6 100644 --- a/Detectors/Upgrades/IT3/README.md +++ b/Detectors/Upgrades/IT3/README.md @@ -18,4 +18,4 @@ o2-sim -m PIPE IT3 [...] ``` ## Reconstruction -Currently, the reconstruction is driven the `macro/run_trac_alice3.C`, which takes as input the hits generated by simulation. Hits can be smeared upon request (see `kUseSmearing`). \ No newline at end of file +Currently, the reconstruction is driven the `Detectors/Upgrades/ALICE3/TRK/macros/ALICE3toAO2D.C`, which takes as input the hits generated by simulation. Hits can be smeared upon request (see `kUseSmearing`). \ No newline at end of file diff --git a/macro/CMakeLists.txt b/macro/CMakeLists.txt index 97ffe710215e0..5d9a708299020 100644 --- a/macro/CMakeLists.txt +++ b/macro/CMakeLists.txt @@ -52,7 +52,6 @@ install(FILES CheckDigits_mft.C run_rawdecoding_mft.C run_trac_ca_its.C run_trac_its.C - run_trac_alice3.C CreateBCPattern.C UploadDummyAlignment.C DESTINATION share/macro/) @@ -353,11 +352,6 @@ o2_add_test_root_macro(run_trac_its.C O2::SimulationDataFormat LABELS its) -# FIXME: move to subsystem dir -o2_add_test_root_macro(run_trac_alice3.C - PUBLIC_LINK_LIBRARIES O2::GPUTracking O2::ITSMFTSimulation - LABELS its COMPILE_ONLY) - # # NOTE: commented out until unit testing reenabled FIXME : re-enable or delete ? # diff --git a/macro/run_primary_vertexer_ITS.C b/macro/run_primary_vertexer_ITS.C index f1e2ae9176d32..8dc43d70d4d00 100644 --- a/macro/run_primary_vertexer_ITS.C +++ b/macro/run_primary_vertexer_ITS.C @@ -187,18 +187,19 @@ int run_primary_vertexer_ITS(const GPUDataTypes::DeviceType dtype = GPUDataTypes vertexer.setDebugCentroidsHistograms(); // \debug - total[0] = vertexer.evaluateTask(&o2::its::Vertexer::initialiseVertexer, "Vertexer initialisation", std::cout, eventptr); + auto logger = [](std::string s) { std::cout << s << std::endl; }; + total[0] = vertexer.evaluateTask(&o2::its::Vertexer::initialiseVertexer, "Vertexer initialisation", logger, eventptr); // total[1] = vertexer.evaluateTask(&o2::its::Vertexer::findTrivialMCTracklets, "Trivial Tracklet finding", std::cout); // If enable this, comment out the validateTracklets - total[1] = vertexer.evaluateTask(&o2::its::Vertexer::findTracklets, "Tracklet finding", std::cout); + total[1] = vertexer.evaluateTask(&o2::its::Vertexer::findTracklets, "Tracklet finding", logger); #ifdef _ALLOW_DEBUG_TREES_ITS_ if (useMCcheck) { vertexer.evaluateTask(&o2::its::Vertexer::filterMCTracklets, "MC tracklets filtering", std::cout); } #endif - total[2] = vertexer.evaluateTask(&o2::its::Vertexer::validateTracklets, "Adjacent tracklets validation", std::cout); + total[2] = vertexer.evaluateTask(&o2::its::Vertexer::validateTracklets, "Adjacent tracklets validation", logger); // In case willing to use the histogram-based CPU vertexer // total[3] = vertexer.evaluateTask(&o2::its::Vertexer::findHistVertices, "Vertex finding with histograms", std::cout); - total[3] = vertexer.evaluateTask(&o2::its::Vertexer::findVertices, "Vertex finding", std::cout); + total[3] = vertexer.evaluateTask(&o2::its::Vertexer::findVertices, "Vertex finding", logger); std::vector vertITS = vertexer.exportVertices(); const size_t numVert = vertITS.size(); diff --git a/macro/run_trac_alice3.C b/macro/run_trac_alice3.C deleted file mode 100644 index 8fe3f7af9da37..0000000000000 --- a/macro/run_trac_alice3.C +++ /dev/null @@ -1,263 +0,0 @@ -#if !defined(__CLING__) || defined(__ROOTCLING__) -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include "DataFormatsITSMFT/ROFRecord.h" -#include "ITSMFTSimulation/Hit.h" -#include "ITStracking/Configuration.h" -#include "ITStracking/IOUtils.h" -#include "ITStracking/Tracker.h" -#include "ITStracking/TrackerTraitsCPU.h" -#include "ITStracking/Vertexer.h" -#include "DataFormatsITSMFT/Cluster.h" -#include "DataFormatsITSMFT/TopologyDictionary.h" -#include "DataFormatsITSMFT/CompCluster.h" -#include "DetectorsCommonDataFormats/DetID.h" -#include "SimulationDataFormat/MCTrack.h" -#include "MathUtils/Cartesian.h" -#include "ReconstructionDataFormats/DCA.h" -#include "ReconstructionDataFormats/Vertex.h" -#endif - -using o2::its::MemoryParameters; -using o2::its::TrackingParameters; -using o2::itsmft::Hit; -using std::string; - -constexpr bool kUseSmearing{false}; - -struct particle { - int pdg = 0; - int nLayers = 0; - float pt; - float eta; - float phi; - float recoPt; - float recoEta; - float energyFirst; - float energyLast; - int isReco = 0; // 0 = no, 1 = good, 2 = fake -}; - -float getDetLengthFromEta(const float eta, const float radius) -{ - return 10. * (10. + radius * std::cos(2 * std::atan(std::exp(-eta)))); -} - -void run_trac_alice3(const string hitsFileName = "o2sim_HitsTRK.root") -{ - - TChain mcTree("o2sim"); - mcTree.AddFile("o2sim_Kine.root"); - mcTree.SetBranchStatus("*", 0); //disable all branches - mcTree.SetBranchStatus("MCTrack*", 1); - - std::vector* mcArr = nullptr; - mcTree.SetBranchAddress("MCTrack", &mcArr); - - o2::its::Vertexer vertexer(new o2::its::VertexerTraits()); - TChain itsHits("o2sim"); - - itsHits.AddFile(hitsFileName.data()); - - o2::its::Tracker tracker(new o2::its::TrackerTraitsCPU()); - tracker.setBz(5.f); - tracker.setCorrType(o2::base::PropagatorImpl::MatCorrType::USEMatCorrTGeo); - - std::uint32_t roFrame; - std::vector* hits = nullptr; - itsHits.SetBranchAddress("TRKHit", &hits); - - std::vector trackParams(4); - trackParams[0].NLayers = 10; - trackParams[0].MinTrackLength = 10; - - std::vector LayerRadii = {1.8f, 2.8f, 3.8f, 8.0f, 20.0f, 25.0f, 40.0f, 55.f, 80.0f, 100.f}; - std::vector LayerZ(10); - for (int i{0}; i < 10; ++i) - LayerZ[i] = getDetLengthFromEta(1.44, LayerRadii[i]) + 1.; - std::vector TrackletMaxDeltaZ = {0.1f, 0.1f, 0.3f, 0.3f, 0.3f, 0.3f, 0.5f, 0.5f, 0.5f}; - std::vector CellMaxDCA = {0.05f, 0.04f, 0.05f, 0.2f, 0.4f, 0.5f, 0.5f, 0.5f}; - std::vector CellMaxDeltaZ = {0.2f, 0.4f, 0.5f, 0.6f, 3.0f, 3.0f, 3.0f, 3.0f}; - std::vector NeighbourMaxDeltaCurvature = {0.008f, 0.0025f, 0.003f, 0.0035f, 0.004f, 0.004f, 0.005f}; - std::vector NeighbourMaxDeltaN = {0.002f, 0.0090f, 0.002f, 0.005f, 0.005f, 0.005f, 0.005f}; - - trackParams[0].LayerRadii = LayerRadii; - trackParams[0].LayerZ = LayerZ; - trackParams[0].TrackletMaxDeltaPhi = 0.3; - trackParams[0].CellMaxDeltaPhi = 0.15; - trackParams[0].CellMaxDeltaTanLambda = 0.03; - trackParams[0].TrackletMaxDeltaZ = TrackletMaxDeltaZ; - trackParams[0].CellMaxDCA = CellMaxDCA; - trackParams[0].CellMaxDeltaZ = CellMaxDeltaZ; - trackParams[0].NeighbourMaxDeltaCurvature = NeighbourMaxDeltaCurvature; - trackParams[0].NeighbourMaxDeltaN = NeighbourMaxDeltaN; - - std::vector memParams(4); - std::vector CellsMemoryCoefficients = {2.3208e-08f * 20, 2.104e-08f * 20, 1.6432e-08f * 20, 1.2412e-08f * 20, 1.3543e-08f * 20, 1.5e-08f * 20, 1.6e-08f * 20, 1.7e-08f * 20}; - std::vector TrackletsMemoryCoefficients = {0.0016353f * 1000, 0.0013627f * 1000, 0.000984f * 1000, 0.00078135f * 1000, 0.00057934f * 1000, 0.00052217f * 1000, 0.00052217f * 1000, 0.00052217f * 1000, 0.00052217f * 1000}; - memParams[0].CellsMemoryCoefficients = CellsMemoryCoefficients; - memParams[0].TrackletsMemoryCoefficients = TrackletsMemoryCoefficients; - memParams[0].MemoryOffset = 8000; - - for (int i = 1; i < 4; ++i) { - memParams[i] = memParams[i - 1]; - trackParams[i] = trackParams[i - 1]; - // trackParams[i].MinTrackLength -= 2; - trackParams[i].TrackletMaxDeltaPhi = trackParams[i].TrackletMaxDeltaPhi * 3 > TMath::Pi() ? TMath::Pi() : trackParams[i].TrackletMaxDeltaPhi * 3; - trackParams[i].CellMaxDeltaPhi = trackParams[i].CellMaxDeltaPhi * 3 > TMath::Pi() ? TMath::Pi() : trackParams[i].CellMaxDeltaPhi * 3; - trackParams[i].CellMaxDeltaTanLambda *= 3; - for (auto& val : trackParams[i].TrackletMaxDeltaZ) - val *= 3; - for (auto& val : trackParams[i].CellMaxDCA) - val *= 3; - for (auto& val : trackParams[i].CellMaxDeltaZ) - val *= 3; - for (auto& val : trackParams[i].NeighbourMaxDeltaCurvature) - val *= 3; - for (auto& val : trackParams[i].NeighbourMaxDeltaN) - val *= 3; - } - - tracker.setParameters(memParams, trackParams); - - constexpr int nBins = 100; - constexpr float minPt = 0.01; - constexpr float maxPt = 10; - double newBins[nBins + 1]; - newBins[0] = minPt; - double factor = pow(maxPt / minPt, 1. / nBins); - for (int i = 1; i <= nBins; i++) { - newBins[i] = factor * newBins[i - 1]; - } - - TH1D genH("gen", ";#it{p}_{T} (GeV/#it{c});", nBins, newBins); - TH1D recH("rec", "Efficiency;#it{p}_{T} (GeV/#it{c});", nBins, newBins); - TH1D fakH("fak", "Fake rate;#it{p}_{T} (GeV/#it{c});", nBins, newBins); - TH2D deltaPt("deltapt", ";#it{p}_{T} (GeV/#it{c});#Delta#it{p}_{T} (GeV/#it{c});", nBins, newBins, 100, -0.1, 0.1); - TH2D dcaxy("dcaxy", ";#it{p}_{T} (GeV/#it{c});DCA_{xy} (#mum);", nBins, newBins, 200, -200, 200); - TH2D dcaz("dcaz", ";#it{p}_{T} (GeV/#it{c});DCA_{z} (#mum);", nBins, newBins, 200, -200, 200); - TH2D deltaE("deltaE", ";#it{p}_{T} (GeV/#it{c});#DeltaE (MeV);", nBins, newBins, 200, 0, 100); - - ROOT::Math::XYZPointF pos{0.f, 0.f, 0.f}; - const std::array cov{1.e-4, 0., 0., 1.e-4, 0., 1.e-4}; - o2::dataformats::VertexBase vtx(pos, cov); - o2::dataformats::DCA dca; - - for (int iEvent{0}; iEvent < itsHits.GetEntriesFast(); ++iEvent) { - std::cout << "*************** Event " << iEvent << " ***************" << std::endl; - itsHits.GetEntry(iEvent); - mcTree.GetEvent(iEvent); - o2::its::ROframe event{iEvent, 10}; - - int id{0}; - std::map mapPDG; - for (auto& hit : *hits) { - const int layer{hit.GetDetectorID()}; - float xyz[3]{hit.GetX(), hit.GetY(), hit.GetZ()}; - float r{std::hypot(xyz[0], xyz[1])}; - float phi{std::atan2(-xyz[1], -xyz[0]) + o2::its::constants::math::Pi}; - - if (kUseSmearing) { - phi = gRandom->Gaus(phi, std::asin(0.0005f / r)); - xyz[0] = r * std::cos(phi); - xyz[1] = r * std::sin(phi); - xyz[2] = gRandom->Gaus(xyz[2], 0.0005f); - } - - event.addTrackingFrameInfoToLayer(layer, xyz[0], xyz[1], xyz[2], r, phi, std::array{0.f, xyz[2]}, - std::array{0.0005f * 0.0005f, 0.f, 0.0005f * 0.0005f}); - event.addClusterToLayer(layer, xyz[0], xyz[1], xyz[2], event.getClustersOnLayer(layer).size()); - event.addClusterLabelToLayer(layer, o2::MCCompLabel(hit.GetTrackID(), iEvent, iEvent, false)); - event.addClusterExternalIndexToLayer(layer, id++); - if (mapPDG.find(hit.GetTrackID()) == mapPDG.end()) { - mapPDG[hit.GetTrackID()] = particle(); - mapPDG[hit.GetTrackID()].nLayers |= 1 << layer; - mapPDG[hit.GetTrackID()].pt = std::hypot(hit.GetPx(), hit.GetPy()); - if (hit.GetTrackID() < mcArr->size()) { - auto part = mcArr->at(hit.GetTrackID()); - mapPDG[hit.GetTrackID()].energyFirst = part.GetEnergy(); - mapPDG[hit.GetTrackID()].pdg = part.GetPdgCode(); - mapPDG[hit.GetTrackID()].pt = part.GetPt(); - mapPDG[hit.GetTrackID()].eta = part.GetEta(); - mapPDG[hit.GetTrackID()].phi = part.GetPhi(); - } - } else { - mapPDG[hit.GetTrackID()].nLayers |= 1 << layer; - mapPDG[hit.GetTrackID()].energyLast = hit.GetE(); - } - } - roFrame = iEvent; - - vertexer.clustersToVertices(event); - int nPart10layers{0}; - for (auto& part : mapPDG) { - if (part.second.nLayers == 0x3FF) { - nPart10layers++; - genH.Fill(part.second.pt); - deltaE.Fill(part.second.pt, 1000 * (part.second.energyFirst - part.second.energyLast)); - } - } - tracker.clustersToTracks(event); - auto& tracks = tracker.getTracks(); - auto& tracksLabels = tracker.getTrackLabels(); - std::cout << "**** " << nPart10layers << " particles with 10 layers " << tracks.size() << " tracks" << std::endl; - int good{0}; - for (unsigned int i{0}; i < tracks.size(); ++i) { - auto& lab = tracksLabels[i]; - auto& track = tracks[i]; - int trackID = std::abs(lab.getTrackID()); - if (mapPDG.find(trackID) != mapPDG.end()) { - if (!mapPDG[trackID].pdg) - std::cout << "strange" << std::endl; - mapPDG[trackID].isReco = lab.isFake() ? 2 : 1; - mapPDG[trackID].recoPt = track.getPt(); - mapPDG[trackID].recoEta = track.getEta(); - (lab.isFake() ? fakH : recH).Fill(mapPDG[trackID].pt); - deltaPt.Fill(mapPDG[trackID].pt, (mapPDG[trackID].pt - mapPDG[trackID].recoPt) / mapPDG[trackID].pt); - if (!track.propagateToDCA(vtx, tracker.getBz(), &dca)) { - std::cout << "Propagation failed." << std::endl; - } else { - dcaxy.Fill(mapPDG[trackID].pt, dca.getY() * 1.e4); - dcaz.Fill(mapPDG[trackID].pt, dca.getZ() * 1.e4); - } - } - } - } - - TH1D dcaxy_res(recH); - dcaxy_res.Reset(); - dcaxy_res.SetNameTitle("dcaxy_res", ";#it{p}_{T} (GeV/#it{c});#sigma(DCA_{xy}) (#mum);"); - TH1D dcaz_res(recH); - dcaz_res.Reset(); - dcaz_res.SetNameTitle("dcaz_res", ";#it{p}_{T} (GeV/#it{c});#sigma(DCA_{z}) (#mum);"); - for (int i = 1; i <= nBins; ++i) { - auto proj = dcaxy.ProjectionY(Form("xy%i", i), i, i); - dcaxy_res.SetBinContent(i, proj->GetStdDev()); - dcaxy_res.SetBinError(i, proj->GetStdDevError()); - proj = dcaz.ProjectionY(Form("z%i", i), i, i); - dcaz_res.SetBinContent(i, proj->GetStdDev()); - dcaz_res.SetBinError(i, proj->GetStdDevError()); - } - - TFile output("output.root", "recreate"); - recH.Divide(&genH); - fakH.Divide(&genH); - recH.Write(); - fakH.SetLineColor(kRed); - fakH.Write(); - genH.Write(); - deltaPt.Write(); - dcaxy.Write(); - dcaz.Write(); - dcaz_res.Write(); - dcaxy_res.Write(); - deltaE.Write(); -} diff --git a/macro/run_trac_ca_its.C b/macro/run_trac_ca_its.C index 8806e3b52410d..129adba7bfaf5 100644 --- a/macro/run_trac_ca_its.C +++ b/macro/run_trac_ca_its.C @@ -28,8 +28,10 @@ #include "ITSBase/GeometryTGeo.h" +#include "DataFormatsITSMFT/CompCluster.h" #include "ITStracking/ROframe.h" #include "ITStracking/IOUtils.h" +#include "ITStracking/TimeFrame.h" #include "ITStracking/Tracker.h" #include "ITStracking/TrackerTraitsCPU.h" #include "ITStracking/Vertexer.h" @@ -67,13 +69,6 @@ void run_trac_ca_its(bool cosmics = false, gSystem->Load("libO2ITStracking"); - //std::unique_ptr rec(GPUReconstruction::CreateInstance()); - // std::unique_ptr rec(GPUReconstruction::CreateInstance("CUDA", true)); // for GPU with CUDA - // auto* chainITS = rec->AddChain(); - // rec->Init(); - - // o2::its::Tracker tracker(chainITS->GetITSTrackerTraits()); - o2::its::Tracker tracker(new o2::its::TrackerTraitsCPU()); o2::its::ROframe event(0, 7); if (path.back() != '/') { @@ -92,21 +87,6 @@ void run_trac_ca_its(bool cosmics = false, LOG(FATAL) << "Failed to load ma"; } double origD[3] = {0., 0., 0.}; - tracker.setBz(field->getBz(origD)); - - //-------- init lookuptable --------// - if (useLUT) { - auto* lut = o2::base::MatLayerCylSet::loadFromFile(matLUTFile); - o2::base::Propagator::Instance()->setMatLUT(lut); - } else { - tracker.setCorrType(o2::base::PropagatorImpl::MatCorrType::USEMatCorrTGeo); - } - - if (tracker.isMatLUT()) { - LOG(INFO) << "Loaded material LUT from " << matLUTFile; - } else { - LOG(INFO) << "Material LUT " << matLUTFile << " file is absent, only TGeo can be used"; - } // @@ -231,42 +211,42 @@ void run_trac_ca_its(bool cosmics = false, memParams.resize(3); trackParams[0].TrackletMaxDeltaPhi = 0.05f; trackParams[1].TrackletMaxDeltaPhi = 0.1f; - trackParams[2].MinTrackLength = 4; + trackParams[2].MinTrackLength = 5; trackParams[2].TrackletMaxDeltaPhi = 0.3; // --- - // pp tracking params - // 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 = 7; - // 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); - // for (auto& coef : memParams[1].CellsMemoryCoefficients) - // coef *= 40; - // for (auto& coef : memParams[1].TrackletsMemoryCoefficients) - // coef *= 40; + // 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); // --- } - tracker.setParameters(memParams, trackParams); int currentEvent = -1; gsl::span patt(patterns->data(), patterns->size()); auto pattIt = patt.begin(); auto clSpan = gsl::span(cclusters->data(), cclusters->size()); + o2::its::TimeFrame tf; + gsl::span rofspan(*rofs); + tf.loadROFrameData(rofspan, clSpan, pattIt, dict, labels); + pattIt = patt.begin(); + int rofId{0}; for (auto& rof : *rofs) { auto start = std::chrono::high_resolution_clock::now(); @@ -275,10 +255,11 @@ void run_trac_ca_its(bool cosmics = false, vertexer.initialiseVertexer(&event); vertexer.findTracklets(); - // vertexer.filterMCTracklets(); // to use MC check vertexer.validateTracklets(); vertexer.findVertices(); std::vector vertITS = vertexer.exportVertices(); + rofId++; + tf.addPrimaryVertices(vertITS); auto& vtxROF = vertROFvec.emplace_back(rof); // register entry and number of vertices in the vtxROF.setFirstEntry(vertices.size()); // dedicated ROFRecord vtxROF.setNEntries(vertITS.size()); @@ -295,19 +276,36 @@ void run_trac_ca_its(bool cosmics = false, } trackClIdx.clear(); tracksITS.clear(); - tracker.clustersToTracks(event); - auto end = std::chrono::high_resolution_clock::now(); std::chrono::duration diff_t{end - start}; ncls.push_back(event.getTotalClusters()); time.push_back(diff_t.count()); + roFrameCounter++; + } + tf.printVertices(); - tracks.swap(tracker.getTracks()); - if (tracks.size()) { - std::cout << "\t\tFound " << tracks.size() << " tracks" << std::endl; - } - for (auto& trc : tracks) { + 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); + } else { + tracker.setCorrType(o2::base::PropagatorImpl::MatCorrType::USEMatCorrTGeo); + } + + if (tracker.isMatLUT()) { + LOG(INFO) << "Loaded material LUT from " << matLUTFile; + } else { + LOG(INFO) << "Material LUT " << matLUTFile << " file is absent, only TGeo can be used"; + } + + for (int iROF{0}; iROF < tf.getNrof(); ++iROF) { + tracksITS.clear(); + for (auto& trc : tf.getTracks(iROF)) { trc.setFirstClusterEntry(trackClIdx.size()); // before adding tracks, create final cluster indices int ncl = trc.getNumberOfClusters(); for (int ic = 0; ic < ncl; ic++) { @@ -315,15 +313,13 @@ void run_trac_ca_its(bool cosmics = false, } tracksITS.emplace_back(trc); } - - trackLabels = tracker.getTrackLabels(); /// FIXME: assignment ctor is not optimal. + trackLabels = tf.getTracksLabel(iROF); /// FIXME: assignment ctor is not optimal. outTree.Fill(); - roFrameCounter++; } outFile.cd(); outTree.Write(); - outFile.Close(); + // outFile.Close(); TGraph* graph = new TGraph(ncls.size(), ncls.data(), time.data()); graph->SetMarkerStyle(20);