diff --git a/Detectors/Upgrades/CMakeLists.txt b/Detectors/Upgrades/CMakeLists.txt index 291c1bc465329..0123c1b3de901 100644 --- a/Detectors/Upgrades/CMakeLists.txt +++ b/Detectors/Upgrades/CMakeLists.txt @@ -10,5 +10,5 @@ # or submit itself to any jurisdiction. message(STATUS "Building detectors for upgrades") -add_subdirectory(IT3) +add_subdirectory(ITS3) add_subdirectory(ALICE3) \ No newline at end of file diff --git a/Detectors/Upgrades/IT3/macros/test/CMakeLists.txt b/Detectors/Upgrades/IT3/macros/test/CMakeLists.txt deleted file mode 100644 index db0b8a378d221..0000000000000 --- a/Detectors/Upgrades/IT3/macros/test/CMakeLists.txt +++ /dev/null @@ -1,33 +0,0 @@ -# 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. - -# o2_add_test_root_macro(CheckDigitsITS3.C -# PUBLIC_LINK_LIBRARIES O2::ITSBase -# O2::ITS3Base -# O2::ITSMFTBase -# O2::ITSMFTSimulation -# O2::ITS3Simulation -# O2::MathUtils -# O2::SimulationDataFormat -# O2::DetectorsBase -# LABELS its3) - -# o2_add_test_root_macro(CheckClustersITS3.C -# PUBLIC_LINK_LIBRARIES O2::ITSBase -# O2::ITS3Base -# O2::ITSMFTBase -# O2::ITSMFTSimulation -# O2::ITS3Simulation -# O2::ITS3Reconstruction -# O2::MathUtils -# O2::SimulationDataFormat -# O2::DetectorsBase -# LABELS its3) \ No newline at end of file diff --git a/Detectors/Upgrades/IT3/CMakeLists.txt b/Detectors/Upgrades/ITS3/CMakeLists.txt similarity index 87% rename from Detectors/Upgrades/IT3/CMakeLists.txt rename to Detectors/Upgrades/ITS3/CMakeLists.txt index 2d9ff8a9ac78e..0c91f53afb082 100644 --- a/Detectors/Upgrades/IT3/CMakeLists.txt +++ b/Detectors/Upgrades/ITS3/CMakeLists.txt @@ -9,5 +9,8 @@ # granted to it by virtue of its status as an Intergovernmental Organization # or submit itself to any jurisdiction. +add_subdirectory(macros) add_subdirectory(simulation) add_subdirectory(base) +add_subdirectory(workflow) +add_subdirectory(reconstruction) \ No newline at end of file diff --git a/Detectors/Upgrades/IT3/README.md b/Detectors/Upgrades/ITS3/README.md similarity index 99% rename from Detectors/Upgrades/IT3/README.md rename to Detectors/Upgrades/ITS3/README.md index 1a844ff7e8871..1811bdba48086 100644 --- a/Detectors/Upgrades/IT3/README.md +++ b/Detectors/Upgrades/ITS3/README.md @@ -2,7 +2,7 @@ \page refDetectorsUpgradesIT3 UpgradesIT3 /doxy --> -# IT3 +# ITS3 Upgraded version of the ITS that includes upgraded truly-cylindrical inner barrel. # Run the full simulation diff --git a/Detectors/Upgrades/IT3/base/CMakeLists.txt b/Detectors/Upgrades/ITS3/base/CMakeLists.txt similarity index 100% rename from Detectors/Upgrades/IT3/base/CMakeLists.txt rename to Detectors/Upgrades/ITS3/base/CMakeLists.txt diff --git a/Detectors/Upgrades/IT3/base/include/ITS3Base/MisalignmentParameter.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/MisalignmentParameter.h similarity index 100% rename from Detectors/Upgrades/IT3/base/include/ITS3Base/MisalignmentParameter.h rename to Detectors/Upgrades/ITS3/base/include/ITS3Base/MisalignmentParameter.h diff --git a/Detectors/Upgrades/IT3/base/include/ITS3Base/SegmentationSuperAlpide.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationSuperAlpide.h similarity index 92% rename from Detectors/Upgrades/IT3/base/include/ITS3Base/SegmentationSuperAlpide.h rename to Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationSuperAlpide.h index 5605c4f8ec7e9..ab0afae25a1ea 100644 --- a/Detectors/Upgrades/IT3/base/include/ITS3Base/SegmentationSuperAlpide.h +++ b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationSuperAlpide.h @@ -18,6 +18,7 @@ #include #include "MathUtils/Cartesian.h" #include "CommonConstants/MathConstants.h" +#include namespace o2 { @@ -28,19 +29,18 @@ namespace its3 class SegmentationSuperAlpide { public: - SegmentationSuperAlpide(int layer = 0) : mLayer{layer}, - NPixels{NRows * NCols}, - NRows{static_cast(double(Radii[layer]) * double(constants::math::TwoPI) / double(PitchRow) + 1)}, + SegmentationSuperAlpide(int layer = 0) : NPixels{NRows * NCols}, + NRows{static_cast(double(Radii[layer]) * double(constants::math::PI) / double(PitchRow) + 1)}, ActiveMatrixSizeRows{PitchRow * NRows}, SensorSizeRows{ActiveMatrixSizeRows + PassiveEdgeTop + PassiveEdgeReadOut} { + LOGP(info, "rows: {} cols: {} npixels: {}", NRows, NCols, NPixels); + LOGP(info, "SegmentationSuperAlpide: layer {} ActiveMatrixSizeRows: {} ActiveMatrixSizeCols: {}", layer, ActiveMatrixSizeCols, ActiveMatrixSizeRows); } - int mLayer; - static constexpr int NLayers = 4; + static constexpr std::array Radii = {1.8f, 2.4f, 3.0f, 7.0f, 10.f}; static constexpr float Length = 27.15f; - static constexpr float Radii[NLayers] = {1.8f, 2.4f, 3.0f, 7.0f}; - static constexpr float PitchCol = 29.24e-4; - static constexpr float PitchRow = 26.88e-4; + static constexpr float PitchCol = 20.e-4; + static constexpr float PitchRow = 20.e-4; static constexpr int NCols = Length / PitchCol; int NRows; int NPixels; diff --git a/Detectors/Upgrades/IT3/base/src/ITS3BaseLinkDef.h b/Detectors/Upgrades/ITS3/base/src/ITS3BaseLinkDef.h similarity index 100% rename from Detectors/Upgrades/IT3/base/src/ITS3BaseLinkDef.h rename to Detectors/Upgrades/ITS3/base/src/ITS3BaseLinkDef.h diff --git a/Detectors/Upgrades/IT3/base/src/MisalignmentParameter.cxx b/Detectors/Upgrades/ITS3/base/src/MisalignmentParameter.cxx similarity index 100% rename from Detectors/Upgrades/IT3/base/src/MisalignmentParameter.cxx rename to Detectors/Upgrades/ITS3/base/src/MisalignmentParameter.cxx diff --git a/Detectors/Upgrades/IT3/base/src/SegmentationSuperAlpide.cxx b/Detectors/Upgrades/ITS3/base/src/SegmentationSuperAlpide.cxx similarity index 95% rename from Detectors/Upgrades/IT3/base/src/SegmentationSuperAlpide.cxx rename to Detectors/Upgrades/ITS3/base/src/SegmentationSuperAlpide.cxx index 10ce7f9a98fb4..ec5e656e36480 100644 --- a/Detectors/Upgrades/IT3/base/src/SegmentationSuperAlpide.cxx +++ b/Detectors/Upgrades/ITS3/base/src/SegmentationSuperAlpide.cxx @@ -21,7 +21,6 @@ using namespace o2::its3; void SegmentationSuperAlpide::print() { - printf("ITS3 segmentation for layer: %d \n", mLayer); printf("Pixel size: %.2f (along %d rows) %.2f (along %d columns) microns\n", PitchRow * 1e4, NRows, PitchCol * 1e4, NCols); printf("Passive edges: bottom: %.2f, top: %.2f, left/right: %.2f microns\n", PassiveEdgeReadOut * 1e4, PassiveEdgeTop * 1e4, PassiveEdgeSide * 1e4); diff --git a/Detectors/Upgrades/IT3/macros/CMakeLists.txt b/Detectors/Upgrades/ITS3/macros/CMakeLists.txt similarity index 100% rename from Detectors/Upgrades/IT3/macros/CMakeLists.txt rename to Detectors/Upgrades/ITS3/macros/CMakeLists.txt diff --git a/Detectors/Upgrades/ITS3/macros/test/CMakeLists.txt b/Detectors/Upgrades/ITS3/macros/test/CMakeLists.txt new file mode 100644 index 0000000000000..4bc7f426bc9cd --- /dev/null +++ b/Detectors/Upgrades/ITS3/macros/test/CMakeLists.txt @@ -0,0 +1,45 @@ +# 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. + +o2_add_test_root_macro(CheckDigitsITS3.C + PUBLIC_LINK_LIBRARIES O2::ITSBase + O2::ITS3Base + O2::ITSMFTBase + O2::ITSMFTSimulation + O2::ITS3Simulation + O2::MathUtils + O2::SimulationDataFormat + O2::DetectorsBase + LABELS its3) + +o2_add_test_root_macro(CheckSquasherITS3.C + PUBLIC_LINK_LIBRARIES O2::ITSBase + O2::ITS3Base + O2::ITSMFTBase + O2::ITSMFTSimulation + O2::ITS3Simulation + O2::ITS3Reconstruction + O2::MathUtils + O2::SimulationDataFormat + O2::DetectorsBase + LABELS its3) + +o2_add_test_root_macro(CheckClustersITS3.C + PUBLIC_LINK_LIBRARIES O2::ITSBase + O2::ITS3Base + O2::ITSMFTBase + O2::ITSMFTSimulation + O2::ITS3Simulation + O2::ITS3Reconstruction + O2::MathUtils + O2::SimulationDataFormat + O2::DetectorsBase + LABELS its3) \ No newline at end of file diff --git a/Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.C b/Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.C new file mode 100644 index 0000000000000..539e5fb055b6f --- /dev/null +++ b/Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.C @@ -0,0 +1,277 @@ +/// \file CheckClusters.C +/// \brief Simple macro to check ITSU clusters + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#define ENABLE_UPGRADES +#include +#include +#include +#include +#include +#include +#include +#include + +#include "DetectorsCommonDataFormats/DetID.h" +#include "ITSMFTBase/SegmentationAlpide.h" +#include "ITS3Base/SegmentationSuperAlpide.h" +#include "ITSBase/GeometryTGeo.h" +#include "DataFormatsITS3/CompCluster.h" +#include "DataFormatsITSMFT/TopologyDictionary.h" +#include "ITS3Reconstruction/TopologyDictionary.h" +#include "ITSMFTSimulation/Hit.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "MathUtils/Cartesian.h" +#include "MathUtils/Utils.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "DetectorsCommonDataFormats/DetectorNameConf.h" +#endif + +void CheckClustersITS3(int nITS3layers = 3, std::string clusfile = "o2clus_it3.root", std::string hitfile = "o2sim_HitsIT3.root", + std::string inputGeom = "", std::string dictfile = "", bool batch = true) +{ + gROOT->SetBatch(batch); + + // we assume that we have 2 chips per layer + const int nChipsPerLayer = 2; + + const int QEDSourceID = 99; // Clusters from this MC source correspond to QED electrons + + using namespace o2::base; + using namespace o2::its; + + using o2::its3::SegmentationSuperAlpide; + using Segmentation = o2::itsmft::SegmentationAlpide; + using o2::its3::CompClusterExt; + using o2::itsmft::Hit; + using ROFRec = o2::itsmft::ROFRecord; + using MC2ROF = o2::itsmft::MC2ROFRecord; + using HitVec = std::vector; + using MC2HITS_map = std::unordered_map; // maps (track_ID<<16 + chip_ID) to entry in the hit vector + + std::vector hitVecPool; + std::vector mc2hitVec; + + TFile fout("CheckClusters.root", "recreate"); + TNtuple nt("ntc", "cluster ntuple", "ev:lab:hlx:hlz:tx:tz:cgx:cgy:cgz:dx:dy:dz:ex:ez:patid:rof:npx:id"); + + // Geometry + o2::base::GeometryManager::loadGeometry(inputGeom); + auto gman = o2::its::GeometryTGeo::Instance(); + gman->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::T2GRot, + o2::math_utils::TransformType::L2G)); // request cached transforms + + // Hits + TFile fileH(hitfile.data()); + TTree* hitTree = (TTree*)fileH.Get("o2sim"); + std::vector* hitArray = nullptr; + hitTree->SetBranchAddress("IT3Hit", &hitArray); + mc2hitVec.resize(hitTree->GetEntries()); + hitVecPool.resize(hitTree->GetEntries(), nullptr); + + // Clusters + TFile fileC(clusfile.data()); + TTree* clusTree = (TTree*)fileC.Get("o2sim"); + std::vector* clusArr = nullptr; + clusTree->SetBranchAddress("IT3ClusterComp", &clusArr); + std::vector* patternsPtr = nullptr; + auto pattBranch = clusTree->GetBranch("IT3ClusterPatt"); + if (pattBranch) { + pattBranch->SetAddress(&patternsPtr); + } + if (dictfile.empty()) { + dictfile = o2::base::DetectorNameConf::getAlpideClusterDictionaryFileName(o2::detectors::DetID::IT3, "", ".bin"); + } + o2::itsmft::TopologyDictionary dict2; + std::ifstream file(dictfile.c_str()); + if (file.good()) { + LOG(info) << "Running with dictionary: " << dictfile.c_str(); + dict2.readFromFile(dictfile); + } else { + LOG(info) << "Running without dictionary !"; + } + o2::its3::TopologyDictionary dict = dict2; + + // ROFrecords + std::vector rofRecVec, *rofRecVecP = &rofRecVec; + clusTree->SetBranchAddress("IT3ClustersROF", &rofRecVecP); + + // Cluster MC labels + o2::dataformats::MCTruthContainer* clusLabArr = nullptr; + std::vector mc2rofVec, *mc2rofVecP = &mc2rofVec; + if (hitTree && clusTree->GetBranch("IT3ClusterMCTruth")) { + clusTree->SetBranchAddress("IT3ClusterMCTruth", &clusLabArr); + clusTree->SetBranchAddress("IT3ClustersMC2ROF", &mc2rofVecP); + } + + clusTree->GetEntry(0); + int nROFRec = (int)rofRecVec.size(); + std::vector mcEvMin(nROFRec, hitTree->GetEntries()); + std::vector mcEvMax(nROFRec, -1); + + // >> build min and max MC events used by each ROF + for (int imc = mc2rofVec.size(); imc--;) { + const auto& mc2rof = mc2rofVec[imc]; + printf("MCRecord: "); + mc2rof.print(); + if (mc2rof.rofRecordID < 0) { + continue; // this MC event did not contribute to any ROF + } + for (int irfd = mc2rof.maxROF - mc2rof.minROF + 1; irfd--;) { + int irof = mc2rof.rofRecordID + irfd; + if (irof >= nROFRec) { + LOG(error) << "ROF=" << irof << " from MC2ROF record is >= N ROFs=" << nROFRec; + } + if (mcEvMin[irof] > imc) { + mcEvMin[irof] = imc; + } + if (mcEvMax[irof] < imc) { + mcEvMax[irof] = imc; + } + } + } + + // << build min and max MC events used by each ROF + auto pattIt = patternsPtr->cbegin(); + for (int irof = 0; irof < nROFRec; irof++) { + const auto& rofRec = rofRecVec[irof]; + + rofRec.print(); + + // >> read and map MC events contributing to this ROF + for (int im = mcEvMin[irof]; im <= mcEvMax[irof]; im++) { + if (!hitVecPool[im]) { + hitTree->SetBranchAddress("IT3Hit", &hitVecPool[im]); + hitTree->GetEntry(im); + auto& mc2hit = mc2hitVec[im]; + const auto* hitArray = hitVecPool[im]; + for (int ih = hitArray->size(); ih--;) { + const auto& hit = (*hitArray)[ih]; + uint64_t key = (uint64_t(hit.GetTrackID()) << 32) + hit.GetDetectorID(); + mc2hit.emplace(key, ih); + } + } + } + // << cache MC events contributing to this ROF + for (int icl = 0; icl < rofRec.getNEntries(); icl++) { + int clEntry = rofRec.getFirstEntry() + icl; // entry of icl-th cluster of this ROF in the vector of clusters + + const auto& cluster = (*clusArr)[clEntry]; + + float errX{0.f}; + float errZ{0.f}; + int npix = 0; + auto pattID = cluster.getPatternID(); + o2::math_utils::Point3D locC; + auto chipID = cluster.getSensorID(); + if (pattID == o2::its3::CompCluster::InvalidPatternID || dict.isGroup(pattID)) { + o2::itsmft::ClusterPattern patt(pattIt); + locC = dict.getClusterCoordinates(chipID, cluster, patt, false); + } else { + locC = dict.getClusterCoordinates(chipID, cluster); + errX = dict.getErrX(pattID); + errZ = dict.getErrZ(pattID); + npix = dict.getNpixels(pattID); + LOGP(info, "I am invalid and I am on chip {}", chipID); + } + // Transformation to the local --> global + auto gloC = gman->getMatrixL2G(chipID)(locC); + if (chipID / nChipsPerLayer < nITS3layers) { + double radius = SegmentationSuperAlpide::Radii[chipID / nChipsPerLayer]; + + bool isTop = !(chipID % nChipsPerLayer); + double phi = locC.X() / radius + (isTop ? -0.5 : +0.5) * (float)TMath::Pi(); + gloC.SetXYZ(-radius * std::cos(phi), -(isTop ? radius * std::sin(phi) - 0.1 / 2 : radius * std::sin(phi) + 0.1 / 2), locC.Z()); + } + + const auto& lab = (clusLabArr->getLabels(clEntry))[0]; + + if (!lab.isValid() || lab.getSourceID() == QEDSourceID) + continue; + + // get MC info + int trID = lab.getTrackID(); + const auto& mc2hit = mc2hitVec[lab.getEventID()]; + const auto* hitArray = hitVecPool[lab.getEventID()]; + uint64_t key = (uint64_t(trID) << 32) + chipID; + auto hitEntry = mc2hit.find(key); + if (hitEntry == mc2hit.end()) { + LOG(error) << "Failed to find MC hit entry for Tr" << trID << " chipID" << chipID; + continue; + } + const auto& hit = (*hitArray)[hitEntry->second]; + // + float dx = 0, dz = 0; + int ievH = lab.getEventID(); + o2::math_utils::Point3D locH, locHsta; + o2::math_utils::Point3D gloH, gloHsta; + + // mean local position of the hit + locH = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); // inverse conversion from global to local + locHsta = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart()); + + if (chipID / nChipsPerLayer < nITS3layers) { + bool isTop = !(chipID % nChipsPerLayer); + float reShiftedY = isTop ? hit.GetPosStart().Y() - 0.1 / 2 : hit.GetPosStart().Y() + 0.1 / 2; + float startPhi{std::atan2(-reShiftedY, -hit.GetPosStart().X()) + (isTop ? (float)TMath::Pi() / 2 : -(float)TMath::Pi() / 2)}; + float reShiftedEndY = isTop ? hit.GetPos().Y() - 0.1 / 2 : hit.GetPos().Y() + 0.1 / 2; + float endPhi{std::atan2(-reShiftedEndY, -hit.GetPos().X()) + (isTop ? (float)TMath::Pi() / 2 : -(float)TMath::Pi() / 2)}; + float deltaY = reShiftedEndY - reShiftedY; + locH.SetXYZ(SegmentationSuperAlpide::Radii[chipID / nChipsPerLayer] * endPhi, 0.f, hit.GetPos().Z()); + locHsta.SetXYZ(SegmentationSuperAlpide::Radii[chipID / nChipsPerLayer] * startPhi, deltaY, hit.GetPosStart().Z()); + } + auto x0 = locHsta.X(), dltx = locH.X() - x0; + auto y0 = locHsta.Y(), dlty = locH.Y() - y0; + auto z0 = locHsta.Z(), dltz = locH.Z() - z0; + + if (chipID / nChipsPerLayer >= nITS3layers) { + auto r = (0.5 * (Segmentation::SensorLayerThickness - Segmentation::SensorLayerThicknessEff) - y0) / dlty; + locH.SetXYZ(x0 + r * dltx, y0 + r * dlty, z0 + r * dltz); + } else { + // not really precise, but okish + locH.SetXYZ(0.5 * (locH.X() + locHsta.X()), 0.5 * (locH.Y() + locHsta.Y()), 0.5 * (locH.Z() + locHsta.Z())); + } + + std::array data = {(float)lab.getEventID(), (float)trID, + locH.X(), locH.Z(), dltx / dlty, dltz / dlty, + gloC.X(), gloC.Y(), gloC.Z(), + locC.X() - locH.X(), locC.Y() - locH.Y(), locC.Z() - locH.Z(), + errX, errZ, (float)pattID, + (float)rofRec.getROFrame(), (float)npix, (float)chipID}; + nt.Fill(data.data()); + } + } + + auto canvCgXCgY = new TCanvas("canvCgXCgY", "", 1600, 1600); + canvCgXCgY->Divide(2, 2); + canvCgXCgY->cd(1); + nt.Draw("cgy:cgx>>h_cgy_vs_cgx_IB(1000, -10, 10, 1000, -10, 10)", "id < 6", "colz"); + canvCgXCgY->cd(2); + nt.Draw("cgy:cgz>>h_cgy_vs_cgz_IB(1000, -15, 15, 1000, -10, 10)", "id < 6", "colz"); + canvCgXCgY->cd(3); + nt.Draw("cgy:cgx>>h_cgy_vs_cgx_OB(1000, -50, 50, 1000, -50, 50)", "id >= 6", "colz"); + canvCgXCgY->cd(4); + nt.Draw("cgy:cgz>>h_cgy_vs_cgz_OB(1000, -100, 100, 1000, -50, 50)", "id >= 6", "colz"); + canvCgXCgY->SaveAs("it3clusters_y_vs_x_vs_z.pdf"); + + auto canvdXdZ = new TCanvas("canvdXdZ", "", 1600, 800); + canvdXdZ->Divide(2, 1); + canvdXdZ->cd(1); + nt.Draw("dx:dz>>h_dx_vs_dz_IB(1000, -0.025, 0.025, 1000, -0.025, 0.025)", "id < 6", "colz"); + canvdXdZ->cd(2); + nt.Draw("dx:dz>>h_dx_vs_dz_OB(1000, -0.025, 0.025, 1000, -0.025, 0.025)", "id >= 6", "colz"); + canvdXdZ->SaveAs("it3clusters_dx_vs_dz.pdf"); + + // auto c1 = new TCanvas("p1", "pullX"); + // c1->cd(); + // c1->SetLogy(); + // nt.Draw("dx/ex", "abs(dx/ex)<10&&patid<10"); + // auto c2 = new TCanvas("p2", "pullZ"); + // c2->cd(); + // c2->SetLogy(); + // nt.Draw("dz/ez", "abs(dz/ez)<10&&patid<10"); + + fout.cd(); + nt.Write(); +} diff --git a/Detectors/Upgrades/IT3/macros/test/CheckClustersITS3.notest b/Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.notest similarity index 69% rename from Detectors/Upgrades/IT3/macros/test/CheckClustersITS3.notest rename to Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.notest index 7c8d55eaa7b51..539e5fb055b6f 100644 --- a/Detectors/Upgrades/IT3/macros/test/CheckClustersITS3.notest +++ b/Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.notest @@ -9,6 +9,8 @@ #include #include #include +#include +#include #include "DetectorsCommonDataFormats/DetID.h" #include "ITSMFTBase/SegmentationAlpide.h" @@ -26,9 +28,14 @@ #include "DetectorsCommonDataFormats/DetectorNameConf.h" #endif -void CheckClustersITS3(std::string clusfile = "o2clus_it3.root", std::string hitfile = "o2sim_HitsIT3.root", - std::string inputGeom = "", std::string dictfile = "") +void CheckClustersITS3(int nITS3layers = 3, std::string clusfile = "o2clus_it3.root", std::string hitfile = "o2sim_HitsIT3.root", + std::string inputGeom = "", std::string dictfile = "", bool batch = true) { + gROOT->SetBatch(batch); + + // we assume that we have 2 chips per layer + const int nChipsPerLayer = 2; + const int QEDSourceID = 99; // Clusters from this MC source correspond to QED electrons using namespace o2::base; @@ -158,7 +165,7 @@ void CheckClustersITS3(std::string clusfile = "o2clus_it3.root", std::string hit auto pattID = cluster.getPatternID(); o2::math_utils::Point3D locC; auto chipID = cluster.getSensorID(); - if (pattID == o2::itsmft::CompCluster::InvalidPatternID || dict.isGroup(pattID)) { + if (pattID == o2::its3::CompCluster::InvalidPatternID || dict.isGroup(pattID)) { o2::itsmft::ClusterPattern patt(pattIt); locC = dict.getClusterCoordinates(chipID, cluster, patt, false); } else { @@ -166,13 +173,16 @@ void CheckClustersITS3(std::string clusfile = "o2clus_it3.root", std::string hit errX = dict.getErrX(pattID); errZ = dict.getErrZ(pattID); npix = dict.getNpixels(pattID); + LOGP(info, "I am invalid and I am on chip {}", chipID); } // Transformation to the local --> global - auto gloC = gman->getMatrixL2G(chipID) * locC; - if (chipID < o2::its3::SegmentationSuperAlpide::NLayers) { - double radius = SegmentationSuperAlpide::Radii[chipID]; - double phi = locC.X() / radius; - gloC.SetXYZ(radius * std::cos(phi), radius * std::sin(phi), locC.Z()); + auto gloC = gman->getMatrixL2G(chipID)(locC); + if (chipID / nChipsPerLayer < nITS3layers) { + double radius = SegmentationSuperAlpide::Radii[chipID / nChipsPerLayer]; + + bool isTop = !(chipID % nChipsPerLayer); + double phi = locC.X() / radius + (isTop ? -0.5 : +0.5) * (float)TMath::Pi(); + gloC.SetXYZ(-radius * std::cos(phi), -(isTop ? radius * std::sin(phi) - 0.1 / 2 : radius * std::sin(phi) + 0.1 / 2), locC.Z()); } const auto& lab = (clusLabArr->getLabels(clEntry))[0]; @@ -195,24 +205,34 @@ void CheckClustersITS3(std::string clusfile = "o2clus_it3.root", std::string hit float dx = 0, dz = 0; int ievH = lab.getEventID(); o2::math_utils::Point3D locH, locHsta; + o2::math_utils::Point3D gloH, gloHsta; // mean local position of the hit locH = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); // inverse conversion from global to local locHsta = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart()); - if (chipID < 4) { - float startPhi{std::atan2(-hit.GetPosStart().Y(), -hit.GetPosStart().X())}; - float endPhi{std::atan2(-hit.GetPos().Y(), -hit.GetPos().X())}; - locH.SetXYZ(SegmentationSuperAlpide::Radii[chipID] * endPhi, 0.f, hit.GetPos().Z()); - locHsta.SetXYZ(SegmentationSuperAlpide::Radii[chipID] * startPhi, 0.f, hit.GetPosStart().Z()); + if (chipID / nChipsPerLayer < nITS3layers) { + bool isTop = !(chipID % nChipsPerLayer); + float reShiftedY = isTop ? hit.GetPosStart().Y() - 0.1 / 2 : hit.GetPosStart().Y() + 0.1 / 2; + float startPhi{std::atan2(-reShiftedY, -hit.GetPosStart().X()) + (isTop ? (float)TMath::Pi() / 2 : -(float)TMath::Pi() / 2)}; + float reShiftedEndY = isTop ? hit.GetPos().Y() - 0.1 / 2 : hit.GetPos().Y() + 0.1 / 2; + float endPhi{std::atan2(-reShiftedEndY, -hit.GetPos().X()) + (isTop ? (float)TMath::Pi() / 2 : -(float)TMath::Pi() / 2)}; + float deltaY = reShiftedEndY - reShiftedY; + locH.SetXYZ(SegmentationSuperAlpide::Radii[chipID / nChipsPerLayer] * endPhi, 0.f, hit.GetPos().Z()); + locHsta.SetXYZ(SegmentationSuperAlpide::Radii[chipID / nChipsPerLayer] * startPhi, deltaY, hit.GetPosStart().Z()); } - auto x0 = locHsta.X(), dltx = locH.X() - x0; auto y0 = locHsta.Y(), dlty = locH.Y() - y0; auto z0 = locHsta.Z(), dltz = locH.Z() - z0; - auto r = (0.5 * (Segmentation::SensorLayerThickness - Segmentation::SensorLayerThicknessEff) - y0) / dlty; - locH.SetXYZ(x0 + r * dltx, y0 + r * dlty, z0 + r * dltz); - //locH.SetXYZ(0.5 * (locH.X() + locHsta.X()), 0.5 * (locH.Y() + locHsta.Y()), 0.5 * (locH.Z() + locHsta.Z())); + + if (chipID / nChipsPerLayer >= nITS3layers) { + auto r = (0.5 * (Segmentation::SensorLayerThickness - Segmentation::SensorLayerThicknessEff) - y0) / dlty; + locH.SetXYZ(x0 + r * dltx, y0 + r * dlty, z0 + r * dltz); + } else { + // not really precise, but okish + locH.SetXYZ(0.5 * (locH.X() + locHsta.X()), 0.5 * (locH.Y() + locHsta.Y()), 0.5 * (locH.Z() + locHsta.Z())); + } + std::array data = {(float)lab.getEventID(), (float)trID, locH.X(), locH.Z(), dltx / dlty, dltz / dlty, gloC.X(), gloC.Y(), gloC.Z(), @@ -223,30 +243,34 @@ void CheckClustersITS3(std::string clusfile = "o2clus_it3.root", std::string hit } } - new TCanvas; - nt.Draw("cgy:cgx"); - new TCanvas; - nt.Draw("dz:dx", "abs(dz)<0.01 && abs(dx)<0.01"); - new TCanvas; - nt.Draw("dz:tz", "abs(dz)<0.005 && abs(tz)<2"); - - auto c1 = new TCanvas("p1", "pullX"); - c1->cd(); - c1->SetLogy(); - nt.Draw("dx/ex", "abs(dx/ex)<10&&patid<10"); - auto c2 = new TCanvas("p2", "pullZ"); - c2->cd(); - c2->SetLogy(); - nt.Draw("dz/ez", "abs(dz/ez)<10&&patid<10"); - - auto d1 = new TCanvas("d1", "deltaX"); - d1->cd(); - d1->SetLogy(); - nt.Draw("dx", "abs(dx)<5"); - auto d2 = new TCanvas("d2", "deltaZ"); - d2->cd(); - d2->SetLogy(); - nt.Draw("dz", "abs(dz)<5"); + auto canvCgXCgY = new TCanvas("canvCgXCgY", "", 1600, 1600); + canvCgXCgY->Divide(2, 2); + canvCgXCgY->cd(1); + nt.Draw("cgy:cgx>>h_cgy_vs_cgx_IB(1000, -10, 10, 1000, -10, 10)", "id < 6", "colz"); + canvCgXCgY->cd(2); + nt.Draw("cgy:cgz>>h_cgy_vs_cgz_IB(1000, -15, 15, 1000, -10, 10)", "id < 6", "colz"); + canvCgXCgY->cd(3); + nt.Draw("cgy:cgx>>h_cgy_vs_cgx_OB(1000, -50, 50, 1000, -50, 50)", "id >= 6", "colz"); + canvCgXCgY->cd(4); + nt.Draw("cgy:cgz>>h_cgy_vs_cgz_OB(1000, -100, 100, 1000, -50, 50)", "id >= 6", "colz"); + canvCgXCgY->SaveAs("it3clusters_y_vs_x_vs_z.pdf"); + + auto canvdXdZ = new TCanvas("canvdXdZ", "", 1600, 800); + canvdXdZ->Divide(2, 1); + canvdXdZ->cd(1); + nt.Draw("dx:dz>>h_dx_vs_dz_IB(1000, -0.025, 0.025, 1000, -0.025, 0.025)", "id < 6", "colz"); + canvdXdZ->cd(2); + nt.Draw("dx:dz>>h_dx_vs_dz_OB(1000, -0.025, 0.025, 1000, -0.025, 0.025)", "id >= 6", "colz"); + canvdXdZ->SaveAs("it3clusters_dx_vs_dz.pdf"); + + // auto c1 = new TCanvas("p1", "pullX"); + // c1->cd(); + // c1->SetLogy(); + // nt.Draw("dx/ex", "abs(dx/ex)<10&&patid<10"); + // auto c2 = new TCanvas("p2", "pullZ"); + // c2->cd(); + // c2->SetLogy(); + // nt.Draw("dz/ez", "abs(dz/ez)<10&&patid<10"); fout.cd(); nt.Write(); diff --git a/Detectors/Upgrades/IT3/macros/test/CheckDigitsITS3.notest b/Detectors/Upgrades/ITS3/macros/test/CheckDigitsITS3.C similarity index 65% rename from Detectors/Upgrades/IT3/macros/test/CheckDigitsITS3.notest rename to Detectors/Upgrades/ITS3/macros/test/CheckDigitsITS3.C index bac555eec6e2e..c41c084c34f31 100644 --- a/Detectors/Upgrades/IT3/macros/test/CheckDigitsITS3.notest +++ b/Detectors/Upgrades/ITS3/macros/test/CheckDigitsITS3.C @@ -9,9 +9,10 @@ #include #include #include +#include #include -#include "ITS3Base/GeometryTGeo.h" +#include "ITSBase/GeometryTGeo.h" #include "DataFormatsITSMFT/Digit.h" #include "ITS3Base/SegmentationSuperAlpide.h" #include "ITSMFTBase/SegmentationAlpide.h" @@ -26,8 +27,9 @@ #endif -void CheckDigitsITS3(std::string digifile = "it3digits.root", std::string hitfile = "o2sim_HitsIT3.root", std::string inputGeom = "", std::string paramfile = "o2sim_par.root") +void CheckDigitsITS3(int nITS3layers = 3, std::string digifile = "it3digits.root", std::string hitfile = "o2sim_HitsIT3.root", std::string inputGeom = "", std::string paramfile = "o2sim_par.root", bool batch = true) { + gROOT->SetBatch(batch); using namespace o2::base; @@ -38,16 +40,22 @@ void CheckDigitsITS3(std::string digifile = "it3digits.root", std::string hitfil using o2::itsmft::SegmentationAlpide; TFile* f = TFile::Open("CheckDigits.root", "recreate"); - TNtuple* nt = new TNtuple("ntd", "digit ntuple", "id:x:y:z:rowD:colD:rowH:colH:xlH:zlH:xlcH:zlcH:dx:dz"); // Geometry o2::base::GeometryManager::loadGeometry(inputGeom); - auto* gman = o2::its3::GeometryTGeo::Instance(); + auto* gman = o2::its::GeometryTGeo::Instance(); gman->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); - SegmentationSuperAlpide segs[4]{SegmentationSuperAlpide(0), SegmentationSuperAlpide(1), SegmentationSuperAlpide(2), SegmentationSuperAlpide(3)}; - SegmentationSuperAlpide& seg = segs[0]; + // we assume that we have 2 chips per layer + const int nChipsPerLayer = 2; + + std::vector segs{}; + for (int iLayer{0}; iLayer < nITS3layers; ++iLayer) { + for (int iChip{0}; iChip < nChipsPerLayer; ++iChip) { + segs.push_back(SegmentationSuperAlpide(iLayer)); + } + } // Hits TFile* hitFile = TFile::Open(hitfile.data()); @@ -153,7 +161,7 @@ void CheckDigitsITS3(std::string digifile = "it3digits.root", std::string hitfil Int_t chipID = (*digArr)[iDigit].getChipIndex(); - if (chipID < 4) { + if (chipID / nChipsPerLayer < nITS3layers) { segs[chipID].detectorToLocal(ix, iz, x, z); } else { SegmentationAlpide::detectorToLocal(ix, iz, x, z); @@ -170,15 +178,15 @@ void CheckDigitsITS3(std::string digifile = "it3digits.root", std::string hitfil nDigitRead++; auto gloD = gman->getMatrixL2G(chipID)(locD); // convert to global - if (chipID < 4) { + if (chipID / nChipsPerLayer < nITS3layers) { // - // invert - //xyzLocS = {SegmentationSuperAlpide::Radii[detID] * startPhi, 0.f, startPos.Z()}; - //xyzLocE = {SegmentationSuperAlpide::Radii[detID] * endPhi, 0.f, endPos.Z()}; + // invert transformation in Digitiser.cxx // - double radius = SegmentationSuperAlpide::Radii[chipID]; - double phi = locD.X() / radius; - gloD.SetXYZ(radius * std::cos(phi), radius * std::sin(phi), locD.Z()); + double radius = SegmentationSuperAlpide::Radii[chipID / nChipsPerLayer]; + + bool isTop = !(chipID % nChipsPerLayer); + double phi = locD.X() / radius + (isTop ? -0.5 : +0.5) * (float)TMath::Pi(); + gloD.SetXYZ(-radius * std::cos(phi), -(isTop ? radius * std::sin(phi) - 0.1 / 2 : radius * std::sin(phi) + 0.1 / 2), locD.Z()); } float dx = 0., dz = 0.; @@ -194,17 +202,22 @@ void CheckDigitsITS3(std::string digifile = "it3digits.root", std::string hitfil continue; } + ////// HITS Hit& hit = (*hitArray[lab.getEventID()])[hitEntry->second]; // FIXME FOR INNER BARREL auto locH = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); // inverse conversion from global to local auto locHsta = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart()); - if (chipID < 4) { - float startPhi{std::atan2(-hit.GetPosStart().Y(), -hit.GetPosStart().X())}; - float endPhi{std::atan2(-hit.GetPos().Y(), -hit.GetPos().X())}; - locH.SetXYZ(SegmentationSuperAlpide::Radii[chipID] * endPhi, 0.f, hit.GetPos().Z()); - locHsta.SetXYZ(SegmentationSuperAlpide::Radii[chipID] * startPhi, 0.f, hit.GetPosStart().Z()); + if (chipID / nChipsPerLayer < nITS3layers) { + bool isTop = !(chipID % nChipsPerLayer); + float reShiftedY = isTop ? hit.GetPosStart().Y() - 0.1 / 2 : hit.GetPosStart().Y() + 0.1 / 2; + float startPhi{std::atan2(-reShiftedY, -hit.GetPosStart().X()) + (isTop ? (float)TMath::Pi() / 2 : -(float)TMath::Pi() / 2)}; + float reShiftedEndY = isTop ? hit.GetPos().Y() - 0.1 / 2 : hit.GetPos().Y() + 0.1 / 2; + float endPhi{std::atan2(-reShiftedEndY, -hit.GetPos().X()) + (isTop ? (float)TMath::Pi() / 2 : -(float)TMath::Pi() / 2)}; + float deltaY = reShiftedEndY - reShiftedY; + locH.SetXYZ(SegmentationSuperAlpide::Radii[chipID / nChipsPerLayer] * endPhi, 0.f, hit.GetPos().Z()); + locHsta.SetXYZ(SegmentationSuperAlpide::Radii[chipID / nChipsPerLayer] * startPhi, deltaY, hit.GetPosStart().Z()); } locH.SetXYZ(0.5 * (locH.X() + locHsta.X()), 0.5 * (locH.Y() + locHsta.Y()), 0.5 * (locH.Z() + locHsta.Z())); @@ -212,8 +225,13 @@ void CheckDigitsITS3(std::string digifile = "it3digits.root", std::string hitfil int row, col; float xlc = 0., zlc = 0.; - segs[chipID < 4 ? chipID : 0].localToDetector(locH.X(), locH.Z(), row, col); - segs[chipID < 4 ? chipID : 0].detectorToLocal(row, col, xlc, zlc); + if (chipID / nChipsPerLayer < nITS3layers) { + segs[chipID].localToDetector(locH.X(), locH.Z(), row, col); + segs[chipID].detectorToLocal(row, col, xlc, zlc); + } else { + SegmentationAlpide::localToDetector(locH.X(), locH.Z(), row, col); + SegmentationAlpide::detectorToLocal(row, col, xlc, zlc); + } nt->Fill(chipID, gloD.X(), gloD.Y(), gloD.Z(), ix, iz, row, col, locH.X(), locH.Z(), xlc, zlc, locH.X() - locD.X(), locH.Z() - locD.Z()); @@ -224,10 +242,25 @@ void CheckDigitsITS3(std::string digifile = "it3digits.root", std::string hitfil } // end loop on ROFRecords array - new TCanvas; - nt->Draw("y:x"); - new TCanvas; - nt->Draw("dx:dz", "abs(dx)<0.02 && abs(dz)<0.02"); + auto canvXY = new TCanvas("canvXY", "", 1600, 1600); + canvXY->Divide(2, 2); + canvXY->cd(1); + nt->Draw("y:x>>h_y_vs_x_IB(1000, -10, 10, 1000, -10, 10)", "id < 6", "colz"); + canvXY->cd(2); + nt->Draw("y:z>>h_y_vs_z_IB(1000, -15, 15, 1000, -10, 10)", "id < 6", "colz"); + canvXY->cd(3); + nt->Draw("y:x>>h_y_vs_x_OB(1000, -50, 50, 1000, -50, 50)", "id >= 6", "colz"); + canvXY->cd(4); + nt->Draw("y:z>>h_y_vs_z_OB(1000, -100, 100, 1000, -50, 50)", "id >= 6", "colz"); + canvXY->SaveAs("it3digits_y_vs_x_vs_z.pdf"); + + auto canvdXdZ = new TCanvas("canvdXdZ", "", 1600, 800); + canvdXdZ->Divide(2, 1); + canvdXdZ->cd(1); + nt->Draw("dx:dz>>h_dx_vs_dz_IB(1000, -0.025, 0.025, 1000, -0.025, 0.025)", "id < 6", "colz"); + canvdXdZ->cd(2); + nt->Draw("dx:dz>>h_dx_vs_dz_OB(1000, -0.025, 0.025, 1000, -0.025, 0.025)", "id >= 6", "colz"); + canvdXdZ->SaveAs("it3digits_dx_vs_dz.pdf"); f->Write(); f->Close(); diff --git a/Detectors/Upgrades/ITS3/macros/test/CheckSquasherITS3.C b/Detectors/Upgrades/ITS3/macros/test/CheckSquasherITS3.C new file mode 100644 index 0000000000000..b36596abfab21 --- /dev/null +++ b/Detectors/Upgrades/ITS3/macros/test/CheckSquasherITS3.C @@ -0,0 +1,183 @@ +// 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. + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#endif + +void getClusterPatterns(std::vector& pattVec, std::vector* ITSclus, std::vector* ITSpatt, o2::itsmft::TopologyDictionary& mdict); + +void CheckSquasherITS3(const uint chipId = 0, const uint startingROF = 0, const unsigned int nRofs = 3, const string fname = "o2clus_it3.root") +{ + // TColor::InvertPalette(); + gStyle->SetOptStat(0); + gStyle->SetPalette(kInvertedDarkBodyRadiator); + // Geometry + o2::base::GeometryManager::loadGeometry(""); + auto gman = o2::its::GeometryTGeo::Instance(); + gman->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G)); + // Topology dictionary + auto& cc = o2::ccdb::BasicCCDBManager::instance(); + auto mdict = cc.get("ITS/Calib/ClusterDictionary"); + auto fITSclus = TFile::Open(fname.data(), "r"); + auto treeITSclus = (TTree*)fITSclus->Get("o2sim"); + + std::vector* ITSclus = nullptr; + std::vector* ITSrof = nullptr; + std::vector* ITSpatt = nullptr; + o2::dataformats::MCTruthContainer* clusLabArr = nullptr; + + treeITSclus->SetBranchAddress("IT3ClusterComp", &ITSclus); + treeITSclus->SetBranchAddress("IT3ClustersROF", &ITSrof); + treeITSclus->SetBranchAddress("IT3ClusterPatt", &ITSpatt); + // treeITSclus->SetBranchAddress("ITSClusterMCTruth", &clusLabArr); + + auto clSpan = gsl::span(ITSclus->data(), ITSclus->size()); + std::vector hHitMapsVsFrame(nRofs); + + treeITSclus->GetEvent(0); + LOGP(info, "there are {} rofs in this TF", ITSrof->size()); + + // Get patterns + std::vector pattVec; + getClusterPatterns(pattVec, ITSclus, ITSpatt, *mdict); + + for (unsigned int iR{0}; iR < nRofs; iR++) { + LOGP(info, "Processing rof {}", iR + startingROF); + switch(chipId/2) { + case 0: + { + hHitMapsVsFrame[iR] = new TH2D(Form("chip%i_rof%i", chipId, startingROF + iR), Form("chip %i rof %i; ; ; Counts", chipId, startingROF + iR), 13575, -0.5, 13574.5, 2828, -0.5, 2827.5); + break; + } + case 1: + { + hHitMapsVsFrame[iR] = new TH2D(Form("chip%i_rof%i", chipId, startingROF + iR), Form("chip %i rof %i; ; ; Counts", chipId, startingROF + iR), 13575, -0.5, 13574.5, 3770, -0.5, 3769.5); + break; + } + case 2: + { + hHitMapsVsFrame[iR] = new TH2D(Form("chip%i_rof%i", chipId, startingROF + iR), Form("chip %i rof %i; ; ; Counts", chipId, startingROF + iR), 13575, -0.5, 13574.5, 4713, -0.5, 4712.5); + break; + } + default: + { + hHitMapsVsFrame[iR] = new TH2D(Form("chip%i_rof%i", chipId, startingROF + iR), Form("chip %i rof %i; ; ; Counts", chipId, startingROF + iR), 1024, -0.5, 1023.5, 512, -0.5, 511.5); + break; + } + } + + // work on data + const auto& rof = (*ITSrof)[startingROF + iR]; + auto clustersInFrame = rof.getROFData(*ITSclus); + auto patternsInFrame = rof.getROFData(pattVec); + // auto pattIt = ITSpatt->cbegin(); + + for (unsigned int clusInd{0}; clusInd < clustersInFrame.size(); clusInd++) { + const auto& clus = clustersInFrame[clusInd]; + auto sID = clus.getSensorID(); + + if (sID == chipId) { + clus.print(); + + // auto labels = clusLabArr->getLabels(clusInd); + // extract pattern info + auto col = clus.getCol(); + auto row = clus.getRow(); + + std::cout << patternsInFrame[clusInd]; + + std::cout << std::endl; + int ic = 0, ir = 0; + + auto colSpan = patternsInFrame[clusInd].getColumnSpan(); + auto rowSpan = patternsInFrame[clusInd].getRowSpan(); + auto nBits = rowSpan * colSpan; + + for (int i = 2; i < patternsInFrame[clusInd].getUsedBytes() + 2; i++) { + unsigned char tempChar = patternsInFrame[clusInd].getByte(i); + int s = 128; // 0b10000000 + while (s > 0) { + if ((tempChar & s) != 0) // checking active pixels + { + hHitMapsVsFrame[iR]->Fill(col + ic, row + ir); + } + ic++; + s >>= 1; + if ((ir + 1) * ic == nBits) { + break; + } + if (ic == colSpan) { + ic = 0; + ir++; + } + if ((ir + 1) * ic == nBits) { + break; + } + } + } + } + } + } + auto canvas = new TCanvas(Form("chip%d", chipId), Form("chip%d", chipId), nRofs * 1000, 600); + + canvas->Divide(nRofs, 1); + for (unsigned int i{0}; i < nRofs; ++i) { + canvas->cd(i + 1); + gPad->SetGridx(); + gPad->SetGridy(); + hHitMapsVsFrame[i]->Draw("colz"); + } +} + +void getClusterPatterns(std::vector& pattVec, std::vector* ITSclus, std::vector* ITSpatt, o2::itsmft::TopologyDictionary& mdict) +{ + pattVec.reserve(ITSclus->size()); + auto pattIt = ITSpatt->cbegin(); + + for (unsigned int iClus{0}; iClus < ITSclus->size(); ++iClus) { + auto& clus = (*ITSclus)[iClus]; + + auto pattID = clus.getPatternID(); + int npix; + o2::itsmft::ClusterPattern patt; + + if (pattID == o2::its3::CompCluster::InvalidPatternID || mdict.isGroup(pattID)) { + patt.acquirePattern(pattIt); + npix = patt.getNPixels(); + } else { + + npix = mdict.getNpixels(pattID); + patt = mdict.getPattern(pattID); + } + + pattVec.push_back(patt); + } +} diff --git a/Detectors/Upgrades/IT3/reconstruction/CMakeLists.txt b/Detectors/Upgrades/ITS3/reconstruction/CMakeLists.txt similarity index 97% rename from Detectors/Upgrades/IT3/reconstruction/CMakeLists.txt rename to Detectors/Upgrades/ITS3/reconstruction/CMakeLists.txt index e32106dfe261f..dadad61241003 100644 --- a/Detectors/Upgrades/IT3/reconstruction/CMakeLists.txt +++ b/Detectors/Upgrades/ITS3/reconstruction/CMakeLists.txt @@ -13,7 +13,7 @@ o2_add_library(ITS3Reconstruction TARGETVARNAME targetName SOURCES src/Clusterer.cxx src/TopologyDictionary.cxx - + src/FastMultEst.cxx PUBLIC_LINK_LIBRARIES O2::ITSMFTBase O2::ITSMFTReconstruction O2::ITS3Base diff --git a/Detectors/Upgrades/IT3/reconstruction/include/ITS3Reconstruction/Clusterer.h b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h similarity index 63% rename from Detectors/Upgrades/IT3/reconstruction/include/ITS3Reconstruction/Clusterer.h rename to Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h index 9be735ec39e4b..a4f64c4c63f45 100644 --- a/Detectors/Upgrades/IT3/reconstruction/include/ITS3Reconstruction/Clusterer.h +++ b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h @@ -16,6 +16,10 @@ #define _PERFORM_TIMING_ +// uncomment this to not allow diagonal clusters, e.g. like |* | +// | *| +#define _ALLOW_DIAGONAL_ALPIDE_CLUSTERS_ + #include #include #include @@ -51,23 +55,60 @@ class MCTruthContainer; namespace its3 { -using CompClusCont = std::vector; +using CompClusCont = std::vector; using PatternCont = std::vector; using ROFRecCont = std::vector; -//template // container types (PMR or std::vectors) +// template // container types (PMR or std::vectors) class Clusterer { using PixelReader = o2::itsmft::PixelReader; using PixelData = o2::itsmft::PixelData; using ChipPixelData = o2::itsmft::ChipPixelData; + using CompCluster = o2::its3::CompCluster; + using CompClusterExt = o2::its3::CompClusterExt; using Label = o2::MCCompLabel; using MCTruth = o2::dataformats::MCTruthContainer; using ConstMCTruth = o2::dataformats::ConstMCTruthContainerView; public: static constexpr int MaxLabels = 10; + static constexpr int MaxHugeClusWarn = 5; // max number of warnings for HugeCluster + + struct BBox { + uint16_t chipID = 0xffff; + uint16_t rowMin = 0xffff; + uint16_t colMin = 0xffff; + uint16_t rowMax = 0; + uint16_t colMax = 0; + BBox(uint16_t c) : chipID(c) {} + bool isInside(uint16_t row, uint16_t col) const { return row >= rowMin && row <= rowMax && col >= colMin && col <= colMax; } + auto rowSpan() const { return rowMax - rowMin + 1; } + auto colSpan() const { return colMax - colMin + 1; } + bool isAcceptableSize() const { return colMax - colMin < o2::itsmft::ClusterPattern::MaxColSpan && rowMax - rowMin < o2::itsmft::ClusterPattern::MaxRowSpan; } + void clear() + { + rowMin = colMin = 0xffff; + rowMax = colMax = 0; + } + void adjust(uint16_t row, uint16_t col) + { + if (row < rowMin) { + rowMin = row; + } + if (row > rowMax) { + rowMax = row; + } + if (col < colMin) { + colMin = col; + } + if (col > colMax) { + colMax = col; + } + } + }; + //========================================================= /// methods and transient data used within a thread struct ThreadStat { @@ -81,7 +122,7 @@ class Clusterer }; struct ClustererThread { - + int id = -1; Clusterer* parent = nullptr; // parent clusterer // buffers for entries in preClusterIndices in 2 columns, to avoid boundary checks, we reserve // extra elements in the beginning and the end @@ -132,13 +173,8 @@ class Clusterer curr[row] = lastIndex; // store index of the new precluster in the current column buffer } - void streamCluster(const std::vector& pixbuf, uint16_t rowMin, uint16_t rowSpanW, uint16_t colMin, uint16_t colSpanW, - uint16_t chipID, - CompClusCont* compClusPtr, PatternCont* patternsPtr, - MCTruth* labelsClusPtr, int nlab, bool isHuge = false); - void fetchMCLabels(int digID, const ConstMCTruth* labelsDig, int& nfilled); - void initChip(const ChipPixelData* curChipData, uint32_t first, int chipID); + void initChip(const ChipPixelData* curChipData, uint32_t first); void updateChip(const ChipPixelData* curChipData, uint32_t ip); void finishChip(ChipPixelData* curChipData, CompClusCont* compClus, PatternCont* patterns, const ConstMCTruth* labelsDig, MCTruth* labelsClus); @@ -147,7 +183,6 @@ class Clusterer void process(uint16_t chip, uint16_t nChips, CompClusCont* compClusPtr, PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClPtr, const itsmft::ROFRecord& rofPtr); - ClustererThread(Clusterer* par = nullptr) : parent(par) {} ~ClustererThread() { if (column1) { @@ -157,6 +192,11 @@ class Clusterer delete[] column2; } } + ClustererThread(Clusterer* par = nullptr, int _id = -1) : parent(par), id(_id), curr(column2 + 1), prev(column1 + 1) + { + // std::fill(std::begin(column1), std::end(column1), -1); + // std::fill(std::begin(column2), std::end(column2), -1); + } }; //========================================================= @@ -168,6 +208,10 @@ class Clusterer void process(int nThreads, PixelReader& r, CompClusCont* compClus, PatternCont* patterns, ROFRecCont* vecROFRec, MCTruth* labelsCl = nullptr); + template + static void streamCluster(const std::vector& pixbuf, const std::array* lblBuff, const BBox& bbox, const itsmft::LookUp& pattIdConverter, + VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isHuge = false); + bool isContinuousReadOut() const { return mContinuousReadout; } void setContinuousReadOut(bool v) { mContinuousReadout = v; } @@ -177,6 +221,12 @@ class Clusterer int getMaxRowColDiffToMask() const { return mMaxRowColDiffToMask; } void setMaxRowColDiffToMask(int v) { mMaxRowColDiffToMask = v; } + int getMaxROFDepthToSquash() const { return mSquashingDepth; } + void setMaxROFDepthToSquash(int v) { mSquashingDepth = v; } + + int getMaxBCSeparationToSquash() const { return mMaxBCSeparationToSquash; } + void setMaxBCSeparationToSquash(int n) { mMaxBCSeparationToSquash = n; } + void print() const; void clear(); @@ -188,28 +238,12 @@ class Clusterer ///< load the dictionary of cluster topologies void loadDictionary(const std::string& fileName) { mPattIdConverter.loadDictionary(fileName); } + void setDictionary(const itsmft::TopologyDictionary* dict) { mPattIdConverter.setDictionary(dict); } TStopwatch& getTimer() { return mTimer; } // cannot be const TStopwatch& getTimerMerge() { return mTimerMerge; } // cannot be const private: - ///< recalculate min max row and column of the cluster accounting for the position of pix - static void adjustBoundingBox(uint16_t row, uint16_t col, uint16_t& rMin, uint16_t& rMax, uint16_t& cMin, uint16_t& cMax) - { - if (row < rMin) { - rMin = row; - } - if (row > rMax) { - rMax = row; - } - if (col < cMin) { - cMin = col; - } - if (col > cMax) { - cMax = col; - } - } - void flushClusters(CompClusCont* compClus, MCTruth* labels); // clusterization options @@ -218,6 +252,11 @@ class Clusterer ///< mask continuosly fired pixels in frames separated by less than this amount of BCs (fired from hit in prev. ROF) int mMaxBCSeparationToMask = 6000. / o2::constants::lhc::LHCBunchSpacingNS + 10; int mMaxRowColDiffToMask = 0; ///< provide their difference in col/row is <= than this + int mNHugeClus = 0; ///< number of encountered huge clusters + + ///< Squashing options + int mSquashingDepth = 0; ///< squashing is applied to next N rofs + int mMaxBCSeparationToSquash = 6000. / o2::constants::lhc::LHCBunchSpacingNS + 10; std::vector> mThreads; // buffers for threads std::vector mChips; // currently processed ROF's chips data @@ -230,6 +269,51 @@ class Clusterer TStopwatch mTimerMerge; }; +template +void Clusterer::streamCluster(const std::vector& pixbuf, const std::array* lblBuff, const Clusterer::BBox& bbox, const itsmft::LookUp& pattIdConverter, + VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isHuge) +{ + if (labelsClusPtr && lblBuff) { // MC labels were requested + auto cnt = compClusPtr->size(); + for (int i = nlab; i--;) { + labelsClusPtr->addElement(cnt, (*lblBuff)[i]); + } + } + auto colSpanW = bbox.colSpan(); + auto rowSpanW = bbox.rowSpan(); + // add to compact clusters, which must be always filled + std::array patt{}; + for (const auto& pix : pixbuf) { + uint32_t ir = pix.getRowDirect() - bbox.rowMin, ic = pix.getCol() - bbox.colMin; + int nbits = ir * colSpanW + ic; + patt[nbits >> 3] |= (0x1 << (7 - (nbits % 8))); + } + uint16_t pattID = (isHuge || pattIdConverter.size() == 0) ? CompCluster::InvalidPatternID : pattIdConverter.findGroupID(rowSpanW, colSpanW, patt.data()); + uint16_t row = bbox.rowMin, col = bbox.colMin; + if (pattID == CompCluster::InvalidPatternID || pattIdConverter.isGroup(pattID)) { + if (pattID != CompCluster::InvalidPatternID) { + // For groupped topologies, the reference pixel is the COG pixel + float xCOG = 0., zCOG = 0.; + itsmft::ClusterPattern::getCOG(rowSpanW, colSpanW, patt.data(), xCOG, zCOG); + row += round(xCOG); + col += round(zCOG); + } + if (patternsPtr) { + patternsPtr->emplace_back((unsigned char)rowSpanW); + patternsPtr->emplace_back((unsigned char)colSpanW); + int nBytes = rowSpanW * colSpanW / 8; + if (((rowSpanW * colSpanW) % 8) != 0) { + nBytes++; + } + patternsPtr->insert(patternsPtr->end(), std::begin(patt), std::begin(patt) + nBytes); + } + } + compClusPtr->emplace_back(row, col, pattID, bbox.chipID); + if (row > 1024 || col > 1024) { + LOGP(info, "row {} col {}", row, col); + } +} + } // namespace its3 } // namespace o2 #endif /* ALICEO2_ITS_CLUSTERER_H */ \ No newline at end of file diff --git a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/FastMultEst.h b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/FastMultEst.h new file mode 100644 index 0000000000000..058f473c09710 --- /dev/null +++ b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/FastMultEst.h @@ -0,0 +1,66 @@ +// 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. + +#ifndef ALICEO2_ITS3_FASTMULTEST_ +#define ALICEO2_ITS3_FASTMULTEST_ + +#include "ITSMFTReconstruction/ChipMappingITS.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "DataFormatsITS3/CompCluster.h" +#include +#include "ITSReconstruction/FastMultEstConfig.h" +#include +#include + +namespace o2 +{ +namespace its3 +{ + +struct FastMultEst { + + static constexpr int NLayers = o2::itsmft::ChipMappingITS::NLayers; + + float mult = 0.; /// estimated signal clusters multipliciy at reference (1st?) layer + float noisePerChip = 0.; /// estimated or imposed noise per chip + float cov[3] = {0.}; /// covariance matrix of estimation + float chi2 = 0.; /// chi2 + int nLayersUsed = 0; /// number of layers actually used + uint32_t lastRandomSeed = 0; /// state of the gRandom before + + std::array nClPerLayer{0}; // measured N Cl per layer selectROFs + FastMultEst(); + + static uint32_t getCurrentRandomSeed(); + int selectROFs(const gsl::span rofs, const gsl::span clus, + const gsl::span trig, std::vector& sel); + + void fillNClPerLayer(const gsl::span& clusters); + float process(const std::array ncl) + { + return FastMultEstConfig::Instance().imposeNoisePerChip > 0 ? processNoiseImposed(ncl) : processNoiseFree(ncl); + } + float processNoiseFree(const std::array ncl); + float processNoiseImposed(const std::array ncl); + float process(const gsl::span& clusters) + { + fillNClPerLayer(clusters); + return process(nClPerLayer); + } + static bool sSeedSet; + + ClassDefNV(FastMultEst, 1); +}; + +} // namespace its +} // namespace o2 + +#endif diff --git a/Detectors/Upgrades/IT3/reconstruction/include/ITS3Reconstruction/TopologyDictionary.h b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/TopologyDictionary.h similarity index 100% rename from Detectors/Upgrades/IT3/reconstruction/include/ITS3Reconstruction/TopologyDictionary.h rename to Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/TopologyDictionary.h diff --git a/Detectors/Upgrades/IT3/reconstruction/src/Clusterer.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/Clusterer.cxx similarity index 76% rename from Detectors/Upgrades/IT3/reconstruction/src/Clusterer.cxx rename to Detectors/Upgrades/ITS3/reconstruction/src/Clusterer.cxx index 96822b01c178b..5c580672927a8 100644 --- a/Detectors/Upgrades/IT3/reconstruction/src/Clusterer.cxx +++ b/Detectors/Upgrades/ITS3/reconstruction/src/Clusterer.cxx @@ -23,7 +23,6 @@ #ifdef WITH_OPENMP #include #endif - using namespace o2::itsmft; namespace o2::its3 @@ -39,6 +38,7 @@ void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClu nThreads = 1; } auto autoDecode = reader.getDecodeNextAuto(); + int rofcount{0}; do { if (autoDecode) { reader.setDecodeNextAuto(false); // internally do not autodecode @@ -46,6 +46,9 @@ void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClu break; // on the fly decoding was requested, but there were no data left } } + if (reader.getInteractionRecord().isDummy()) { + continue; // No IR info was found + } // pre-fetch all non-empty chips of current ROF ChipPixelData* curChipData = nullptr; mFiredChipsPtr.clear(); @@ -55,7 +58,7 @@ void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClu nPix += curChipData->getData().size(); } - auto& rof = vecROFRec->emplace_back(reader.getInteractionRecord(), 0, compClus->size(), 0); // create new ROF + auto& rof = vecROFRec->emplace_back(reader.getInteractionRecord(), vecROFRec->size(), compClus->size(), 0); // create new ROF uint16_t nFired = mFiredChipsPtr.size(); if (!nFired) { @@ -67,22 +70,20 @@ void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClu if (nFired < nThreads) { nThreads = nFired; } -#ifdef WITH_OPENMP - omp_set_num_threads(nThreads); -#else +#ifndef WITH_OPENMP nThreads = 1; #endif - uint16_t chipStep = nThreads > 1 ? (nThreads == 2 ? 20 : (nThreads < 5 ? 5 : 1)) : nFired; + uint16_t chipStep = nThreads > 1 ? (nThreads == 2 ? 20 : 10) : nFired; int dynGrp = std::min(4, std::max(1, nThreads / 2)); if (nThreads > mThreads.size()) { int oldSz = mThreads.size(); mThreads.resize(nThreads); for (int i = oldSz; i < nThreads; i++) { - mThreads[i] = std::make_unique(this); + mThreads[i] = std::make_unique(this, i); } } #ifdef WITH_OPENMP -#pragma omp parallel for schedule(dynamic, dynGrp) +#pragma omp parallel for schedule(dynamic, dynGrp) num_threads(nThreads) //>> start of MT region for (uint16_t ic = 0; ic < nFired; ic += chipStep) { auto ith = omp_get_thread_num(); @@ -108,6 +109,7 @@ void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClu size_t nClTot = 0, nPattTot = 0; int chid = 0, thrStatIdx[nThreads]; for (int ith = 0; ith < nThreads; ith++) { + std::sort(mThreads[ith]->stats.begin(), mThreads[ith]->stats.end(), [](const ThreadStat& a, const ThreadStat& b) { return a.firstChip < b.firstChip; }); thrStatIdx[ith] = 0; nClTot += mThreads[ith]->compClusters.size(); nPattTot += mThreads[ith]->patterns.size(); @@ -162,10 +164,9 @@ void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClu void Clusterer::ClustererThread::process(uint16_t chip, uint16_t nChips, CompClusCont* compClusPtr, PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClPtr, const ROFRecord& rofPtr) { - if (stats.empty() || stats.back().firstChip + stats.back().nChips < chip) { // there is a jump, register new block + if (stats.empty() || stats.back().firstChip + stats.back().nChips != chip) { // there is a jump, register new block stats.emplace_back(ThreadStat{chip, 0, uint32_t(compClusPtr->size()), patternsPtr ? uint32_t(patternsPtr->size()) : 0, 0, 0}); } - for (int ic = 0; ic < nChips; ic++) { auto* curChipData = parent->mFiredChipsPtr[chip + ic]; auto chipID = curChipData->getChipID(); @@ -175,6 +176,7 @@ void Clusterer::ClustererThread::process(uint16_t chip, uint16_t nChips, CompClu parent->mMaxRowColDiffToMask ? curChipData->maskFiredInSample(parent->mChipsOld[chipID], parent->mMaxRowColDiffToMask) : curChipData->maskFiredInSample(parent->mChipsOld[chipID]); } } + auto nclus0 = compClusPtr->size(); auto validPixID = curChipData->getFirstUnmasked(); auto npix = curChipData->getData().size(); if (validPixID < npix) { // chip data may have all of its pixels masked! @@ -182,7 +184,7 @@ void Clusterer::ClustererThread::process(uint16_t chip, uint16_t nChips, CompClu if (validPixID == npix) { // special case of a single pixel fired on the chip finishChipSingleHitFast(valp, curChipData, compClusPtr, patternsPtr, labelsDigPtr, labelsClPtr); } else { - initChip(curChipData, valp, chipID); + initChip(curChipData, valp); for (; validPixID < npix; validPixID++) { if (!curChipData->getData()[validPixID].isMasked()) { updateChip(curChipData, validPixID); @@ -205,15 +207,13 @@ void Clusterer::ClustererThread::process(uint16_t chip, uint16_t nChips, CompClu void Clusterer::ClustererThread::finishChip(ChipPixelData* curChipData, CompClusCont* compClusPtr, PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClusPtr) { - auto clustersCount = compClusPtr->size(); const auto& pixData = curChipData->getData(); for (int i1 = 0; i1 < preClusterHeads.size(); ++i1) { auto ci = preClusterIndices[i1]; if (ci < 0) { continue; } - uint16_t rowMax = 0, rowMin = 65535; - uint16_t colMax = 0, colMin = 65535; + BBox bbox(curChipData->getChipID()); int nlab = 0; int next = preClusterHeads[i1]; pixArrBuff.clear(); @@ -221,9 +221,13 @@ void Clusterer::ClustererThread::finishChip(ChipPixelData* curChipData, CompClus const auto& pixEntry = pixels[next]; const auto pix = pixData[pixEntry.second]; pixArrBuff.push_back(pix); // needed for cluster topology - adjustBoundingBox(pix.getRowDirect(), pix.getCol(), rowMin, rowMax, colMin, colMax); - if (labelsClusPtr) { // the MCtruth for this pixel is at curChipData->startID+pixEntry.second - fetchMCLabels(pixEntry.second + curChipData->getStartID(), labelsDigPtr, nlab); + bbox.adjust(pix.getRowDirect(), pix.getCol()); + if (labelsClusPtr) { + if (parent->mSquashingDepth) { // the MCtruth for this pixel is stored in chip data: due to squashing we lose contiguity + fetchMCLabels(curChipData->getOrderedPixId(pixEntry.second), labelsDigPtr, nlab); + } else { // the MCtruth for this pixel is at curChipData->startID+pixEntry.second + fetchMCLabels(pixEntry.second + curChipData->getStartID(), labelsDigPtr, nlab); + } } next = pixEntry.first; } @@ -237,107 +241,60 @@ void Clusterer::ClustererThread::finishChip(ChipPixelData* curChipData, CompClus const auto& pixEntry = pixels[next]; const auto pix = pixData[pixEntry.second]; // PixelData pixArrBuff.push_back(pix); // needed for cluster topology - adjustBoundingBox(pix.getRowDirect(), pix.getCol(), rowMin, rowMax, colMin, colMax); - if (labelsClusPtr) { // the MCtruth for this pixel is at curChipData->startID+pixEntry.second - fetchMCLabels(pixEntry.second + curChipData->getStartID(), labelsDigPtr, nlab); + bbox.adjust(pix.getRowDirect(), pix.getCol()); + if (labelsClusPtr) { + if (parent->mSquashingDepth) { // the MCtruth for this pixel is stored in chip data: due to squashing we lose contiguity + fetchMCLabels(curChipData->getOrderedPixId(pixEntry.second), labelsDigPtr, nlab); + } else { // the MCtruth for this pixel is at curChipData->startID+pixEntry.second + fetchMCLabels(pixEntry.second + curChipData->getStartID(), labelsDigPtr, nlab); + } } next = pixEntry.first; } preClusterIndices[i2] = -1; } - - auto chipID = curChipData->getChipID(); - uint16_t colSpan = (colMax - colMin + 1); - uint16_t rowSpan = (rowMax - rowMin + 1); - if (colSpan <= o2::itsmft::ClusterPattern::MaxColSpan && - rowSpan <= o2::itsmft::ClusterPattern::MaxRowSpan) { - streamCluster(pixArrBuff, rowMin, rowSpan, colMin, colSpan, chipID, - compClusPtr, patternsPtr, labelsClusPtr, nlab); + if (bbox.isAcceptableSize()) { + parent->streamCluster(pixArrBuff, &labelsBuff, bbox, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab); } else { - LOG(alarm) << "Splitting a huge cluster ! ChipID: " << chipID; - - colSpan %= o2::itsmft::ClusterPattern::MaxColSpan; - if (colSpan == 0) { - colSpan = o2::itsmft::ClusterPattern::MaxColSpan; - } - - rowSpan %= o2::itsmft::ClusterPattern::MaxRowSpan; - if (rowSpan == 0) { - rowSpan = o2::itsmft::ClusterPattern::MaxRowSpan; + auto warnLeft = MaxHugeClusWarn - parent->mNHugeClus; + if (warnLeft > 0) { + LOGP(warn, "Splitting a huge cluster: chipID {}, rows {}:{} cols {}:{}{}", bbox.chipID, bbox.rowMin, bbox.rowMax, bbox.colMin, bbox.colMax, + warnLeft == 1 ? " (Further warnings will be muted)" : ""); +#ifdef WITH_OPENMP +#pragma omp critical +#endif + { + parent->mNHugeClus++; + } } - + BBox bboxT(bbox); // truncated box + std::vector pixbuf; do { - uint16_t r = rowMin, rsp = rowSpan; - - do { - // Select a subset of pixels fitting the reduced bounding box - std::vector pixbuf; - auto colMax = colMin + colSpan, rowMax = r + rsp; + bboxT.rowMin = bbox.rowMin; + bboxT.colMax = std::min(bbox.colMax, uint16_t(bboxT.colMin + o2::itsmft::ClusterPattern::MaxColSpan - 1)); + do { // Select a subset of pixels fitting the reduced bounding box + bboxT.rowMax = std::min(bbox.rowMax, uint16_t(bboxT.rowMin + o2::itsmft::ClusterPattern::MaxRowSpan - 1)); for (const auto& pix : pixArrBuff) { - if (pix.getRowDirect() >= r && pix.getRowDirect() < rowMax && - pix.getCol() >= colMin && pix.getCol() < colMax) { + if (bboxT.isInside(pix.getRowDirect(), pix.getCol())) { pixbuf.push_back(pix); } } - // Stream a piece of cluster only if the reduced bounding box is not empty - if (!pixbuf.empty()) { - streamCluster(pixbuf, r, rsp, colMin, colSpan, chipID, - compClusPtr, patternsPtr, labelsClusPtr, nlab, true); + if (!pixbuf.empty()) { // Stream a piece of cluster only if the reduced bounding box is not empty + parent->streamCluster(pixbuf, &labelsBuff, bboxT, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab, true); + pixbuf.clear(); } - r += rsp; - rsp = o2::itsmft::ClusterPattern::MaxRowSpan; - } while (r < rowMax); - - colMin += colSpan; - colSpan = o2::itsmft::ClusterPattern::MaxColSpan; - } while (colMin < colMax); + bboxT.rowMin = bboxT.rowMax + 1; + } while (bboxT.rowMin < bbox.rowMax); + bboxT.colMin = bboxT.colMax + 1; + } while (bboxT.colMin < bbox.colMax); } } } -void Clusterer::ClustererThread::streamCluster(const std::vector& pixbuf, uint16_t rowMin, uint16_t rowSpanW, uint16_t colMin, uint16_t colSpanW, uint16_t chipID, CompClusCont* compClusPtr, PatternCont* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isHuge) -{ - if (labelsClusPtr) { // MC labels were requested - auto cnt = compClusPtr->size(); - for (int i = nlab; i--;) { - labelsClusPtr->addElement(cnt, labelsBuff[i]); - } - } - - // add to compact clusters, which must be always filled - unsigned char patt[ClusterPattern::MaxPatternBytes] = {0}; // RSTODO FIX pattern filling - for (const auto& pix : pixbuf) { - unsigned short ir = pix.getRowDirect() - rowMin, ic = pix.getCol() - colMin; - int nbits = ir * colSpanW + ic; - patt[nbits >> 3] |= (0x1 << (7 - (nbits % 8))); - } - uint16_t pattID = (isHuge || parent->mPattIdConverter.size() == 0) ? CompCluster::InvalidPatternID : parent->mPattIdConverter.findGroupID(rowSpanW, colSpanW, patt); - if (pattID == CompCluster::InvalidPatternID || parent->mPattIdConverter.isGroup(pattID)) { - if (pattID != CompCluster::InvalidPatternID) { - //For groupped topologies, the reference pixel is the COG pixel - float xCOG = 0., zCOG = 0.; - ClusterPattern::getCOG(rowSpanW, colSpanW, patt, xCOG, zCOG); - rowMin += round(xCOG); - colMin += round(zCOG); - } - if (patternsPtr) { - patternsPtr->emplace_back((unsigned char)rowSpanW); - patternsPtr->emplace_back((unsigned char)colSpanW); - int nBytes = rowSpanW * colSpanW / 8; - if (((rowSpanW * colSpanW) % 8) != 0) { - nBytes++; - } - patternsPtr->insert(patternsPtr->end(), std::begin(patt), std::begin(patt) + nBytes); - } - } - compClusPtr->emplace_back(rowMin, colMin, pattID, chipID); -} - //__________________________________________________ void Clusterer::ClustererThread::finishChipSingleHitFast(uint32_t hit, ChipPixelData* curChipData, CompClusCont* compClusPtr, PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClusPtr) { - auto clustersCount = compClusPtr->size(); auto pix = curChipData->getData()[hit]; uint16_t row = pix.getRowDirect(), col = pix.getCol(); @@ -373,12 +330,12 @@ Clusterer::Clusterer() : mPattIdConverter() } //__________________________________________________ -void Clusterer::ClustererThread::initChip(const ChipPixelData* curChipData, uint32_t first, int chip) +void Clusterer::ClustererThread::initChip(const ChipPixelData* curChipData, uint32_t first) { // init chip with the 1st unmasked pixel (entry "from" in the mChipData) size = itsmft::SegmentationAlpide::NRows + 2; - if (chip < SegmentationSuperAlpide::NLayers) { - SegmentationSuperAlpide seg(chip); + if (curChipData->getChipID() < 6) { // TODO Fix for mutable layouts + SegmentationSuperAlpide seg(curChipData->getChipID() / 2); size = seg.NRows + 2; } if (column1) { @@ -437,7 +394,11 @@ void Clusterer::ClustererThread::updateChip(const ChipPixelData* curChipData, ui return; } } else { +#ifdef _ALLOW_DIAGONAL_ALPIDE_CLUSTERS_ int neighbours[]{curr[row - 1], prev[row], prev[row + 1], prev[row - 1]}; +#else + int neighbours[]{curr[row - 1], prev[row]}; +#endif for (auto pci : neighbours) { if (pci < 0) { continue; @@ -499,8 +460,10 @@ void Clusterer::clear() void Clusterer::print() const { // print settings + LOGP(info, "Clusterizer squashes overflow pixels separated by {} BC and <= {} in row/col seeking down to {} neighbour ROFs", mMaxBCSeparationToSquash, mMaxRowColDiffToMask, mSquashingDepth); LOG(info) << "Clusterizer masks overflow pixels separated by < " << mMaxBCSeparationToMask << " BC and <= " << mMaxRowColDiffToMask << " in row/col"; + #ifdef _PERFORM_TIMING_ auto& tmr = const_cast(mTimer); // ugly but this is what root does internally auto& tmrm = const_cast(mTimerMerge); @@ -511,5 +474,4 @@ void Clusterer::print() const #endif } - } // namespace o2::its3 \ No newline at end of file diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/FastMultEst.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/FastMultEst.cxx new file mode 100644 index 0000000000000..34928c019446f --- /dev/null +++ b/Detectors/Upgrades/ITS3/reconstruction/src/FastMultEst.cxx @@ -0,0 +1,207 @@ +// 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. + +#include "ITSReconstruction/FastMultEst.h" +#include "ITSMFTBase/DPLAlpideParam.h" +#include "Framework/Logger.h" +#include +#include +#include + +using namespace o2::its; + +bool FastMultEst::sSeedSet = false; + +///______________________________________________________ +FastMultEst::FastMultEst() +{ + if (!sSeedSet && FastMultEstConfig::Instance().cutRandomFraction > 0.f) { + sSeedSet = true; + if (FastMultEstConfig::Instance().randomSeed > 0) { + gRandom->SetSeed(FastMultEstConfig::Instance().randomSeed); + } else if (FastMultEstConfig::Instance().randomSeed < 0) { + gRandom->SetSeed(std::time(nullptr) % 0xffff); + } + } +} + +///______________________________________________________ +/// find multiplicity for given set of clusters +void FastMultEst::fillNClPerLayer(const gsl::span& clusters) +{ + int lr = FastMultEst::NLayers - 1, nchAcc = o2::itsmft::ChipMappingITS::getNChips() - o2::itsmft::ChipMappingITS::getNChipsPerLr(lr); + std::memset(&nClPerLayer[0], 0, sizeof(int) * FastMultEst::NLayers); + for (int i = clusters.size(); i--;) { // profit from clusters being ordered in chip increasing order + while (clusters[i].getSensorID() < nchAcc) { + assert(lr >= 0); + nchAcc -= o2::itsmft::ChipMappingITS::getNChipsPerLr(--lr); + } + nClPerLayer[lr]++; + } +} + +///______________________________________________________ +/// find multiplicity for given number of clusters per layer +float FastMultEst::processNoiseFree(const std::array ncl) +{ + // we assume that on the used layers the observed number of clusters is defined by the + // the noise ~ nu * Nchips and contribution from the signal tracks Ntr*mAccCorr + const auto& conf = FastMultEstConfig::Instance(); + + float mat[3] = {0}, b[2] = {0}; + nLayersUsed = 0; + for (int il = conf.firstLayer; il <= conf.lastLayer; il++) { + if (ncl[il] > 0) { + int nch = o2::itsmft::ChipMappingITS::getNChipsPerLr(il); + float err2i = 1. / ncl[il]; + float m2n = nch * err2i; + mat[0] += err2i * conf.accCorr[il] * conf.accCorr[il]; + mat[2] += nch * m2n; + mat[1] += conf.accCorr[il] * m2n; // non-diagonal element + b[0] += conf.accCorr[il]; + b[1] += nch; + nLayersUsed++; + } + } + mult = noisePerChip = chi2 = -1; + float det = mat[0] * mat[2] - mat[1] * mat[1]; + if (nLayersUsed < 2 || std::abs(det) < 1e-15) { + return -1; + } + float detI = 1. / det; + mult = detI * (b[0] * mat[2] - b[1] * mat[1]); + noisePerChip = detI * (b[1] * mat[0] - b[0] * mat[1]); + cov[0] = mat[2] * detI; + cov[2] = mat[0] * detI; + cov[1] = -mat[1] * detI; + chi2 = 0.; + for (int il = conf.firstLayer; il <= conf.lastLayer; il++) { + if (ncl[il] > 0) { + int nch = o2::itsmft::ChipMappingITS::getNChipsPerLr(il); + float diff = mult * conf.accCorr[il] + nch * noisePerChip - ncl[il]; + chi2 += diff * diff / ncl[il]; + } + } + chi2 = nLayersUsed > 2 ? chi2 / (nLayersUsed - 2) : 0.; + return mult > 0 ? mult : 0; +} + +///______________________________________________________ +/// find multiplicity for given number of clusters per layer with mean noise imposed +float FastMultEst::processNoiseImposed(const std::array ncl) +{ + // we assume that on the used layers the observed number of clusters is defined by the + // the noise ~ nu * Nchips and contribution from the signal tracks Ntr*conf.accCorr + // + // minimize the form sum_lr (noise_i - mu nchips_i)^2 / (mu nchips_i) + lambda_i * (noise_i + mult*acc_i - ncl_i) + // whith noise_i being estimate of the noise clusters in nchips_i of layer i, mu is the mean noise per chip, + // mult is the number of signal clusters on the ref. (1st) layer and the acc_i is the acceptance of layer i wrt 1st. + // The lambda_i is hust a Lagrange multiplier. + + const auto& conf = FastMultEstConfig::Instance(); + float w2sum = 0., wnsum = 0., wsum = 0.; + nLayersUsed = 0; + for (int il = conf.firstLayer; il <= conf.lastLayer; il++) { + if (ncl[il] > 0) { + float nchInv = 1. / o2::itsmft::ChipMappingITS::getNChipsPerLr(il); + w2sum += conf.accCorr[il] * conf.accCorr[il] * nchInv; + wnsum += ncl[il] * nchInv * conf.accCorr[il]; + wsum += conf.accCorr[il]; + nLayersUsed++; + } + } + mult = 0; + chi2 = -1; + noisePerChip = conf.imposeNoisePerChip; + if (nLayersUsed < 1) { + return -1; + } + auto w2sumI = 1. / w2sum; + mult = (wnsum - noisePerChip * wsum) * w2sumI; + cov[0] = wnsum * w2sumI; + cov[2] = 0.; + cov[1] = 0.; + + chi2 = 0.; + for (int il = conf.firstLayer; il <= conf.lastLayer; il++) { + if (ncl[il] > 0) { + float noise = ncl[il] - mult * conf.accCorr[il], estNoise = o2::itsmft::ChipMappingITS::getNChipsPerLr(il) * noisePerChip; + float diff = noise - estNoise; + chi2 += diff * diff / estNoise; + } + } + chi2 = nLayersUsed > 2 ? chi2 / (nLayersUsed - 2) : 0.; + return mult > 0 ? mult : 0; +} + +int FastMultEst::selectROFs(const gsl::span rofs, const gsl::span clus, + const gsl::span trig, std::vector& sel) +{ + int nrof = rofs.size(), nsel = 0; + const auto& multEstConf = FastMultEstConfig::Instance(); // parameters for mult estimation and cuts + sel.clear(); + sel.resize(nrof, true); // by default select all + lastRandomSeed = gRandom->GetSeed(); + if (multEstConf.isMultCutRequested()) { + for (uint32_t irof = 0; irof < nrof; irof++) { + nsel += sel[irof] = multEstConf.isPassingMultCut(process(rofs[irof].getROFData(clus))); + } + } else { + nsel = nrof; + } + using IdNT = std::pair; + if (multEstConf.cutRandomFraction > 0.) { + int ntrig = trig.size(), currTrig = 0; + if (multEstConf.preferTriggered) { + const auto& alpParams = o2::itsmft::DPLAlpideParam::Instance(); + std::vector nTrigROF; + nTrigROF.reserve(nrof); + for (uint32_t irof = 0; irof < nrof; irof++) { + if (sel[irof]) { + if (nsel && gRandom->Rndm() < multEstConf.cutRandomFraction) { + nsel--; + } + auto irROF = rofs[irof].getBCData(); + while (currTrig < ntrig && trig[currTrig].ir < irROF) { // triggers are sorted, jump to 1st one not less than current ROF + currTrig++; + } + auto& trof = nTrigROF.emplace_back(irof, 0); + irROF += alpParams.roFrameLengthInBC; + while (currTrig < ntrig && trig[currTrig].ir < irROF) { + trof.second++; + currTrig++; + } + } + } + if (nsel > 0) { + sort(nTrigROF.begin(), nTrigROF.end(), [](const IdNT& a, const IdNT& b) { return a.second > b.second; }); // order in number of triggers + auto last = nTrigROF.begin() + nsel; + sort(nTrigROF.begin(), last, [](const IdNT& a, const IdNT& b) { return a.first < b.first; }); // order in ROF ID first nsel ROFs + } + for (int i = nsel; i < int(nTrigROF.size()); i++) { // reject ROFs in the tail + sel[nTrigROF[i].first] = false; + } + } else { // dummy random rejection + for (int irof = 0; irof < nrof; irof++) { + if (sel[irof]) { + float sr = gRandom->Rndm(); + if (gRandom->Rndm() < multEstConf.cutRandomFraction) { + sel[irof] = false; + nsel--; + } + } + } + } + } + LOGP(debug, "NSel = {} of {} rofs Seeds: before {} after {}", nsel, nrof, lastRandomSeed, gRandom->GetSeed()); + + return nsel; +} diff --git a/Detectors/Upgrades/IT3/reconstruction/src/ITS3ReconstructionLinkDef.h b/Detectors/Upgrades/ITS3/reconstruction/src/ITS3ReconstructionLinkDef.h similarity index 100% rename from Detectors/Upgrades/IT3/reconstruction/src/ITS3ReconstructionLinkDef.h rename to Detectors/Upgrades/ITS3/reconstruction/src/ITS3ReconstructionLinkDef.h diff --git a/Detectors/Upgrades/IT3/reconstruction/src/TopologyDictionary.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/TopologyDictionary.cxx similarity index 62% rename from Detectors/Upgrades/IT3/reconstruction/src/TopologyDictionary.cxx rename to Detectors/Upgrades/ITS3/reconstruction/src/TopologyDictionary.cxx index ea632ea39401e..51ae9d2699782 100644 --- a/Detectors/Upgrades/IT3/reconstruction/src/TopologyDictionary.cxx +++ b/Detectors/Upgrades/ITS3/reconstruction/src/TopologyDictionary.cxx @@ -13,6 +13,7 @@ #include "ITS3Reconstruction/TopologyDictionary.h" #include "ITS3Base/SegmentationSuperAlpide.h" +#include "ITSMFTBase/SegmentationAlpide.h" namespace o2::its3 { @@ -20,10 +21,14 @@ namespace o2::its3 math_utils::Point3D TopologyDictionary::getClusterCoordinates(int detID, const its3::CompCluster& cl) const { - static SegmentationSuperAlpide segmentations[SegmentationSuperAlpide::NLayers]{SegmentationSuperAlpide(0), SegmentationSuperAlpide(1), SegmentationSuperAlpide(2), SegmentationSuperAlpide(3)}; + static SegmentationSuperAlpide segmentations[6]{SegmentationSuperAlpide(0), + SegmentationSuperAlpide(0), + SegmentationSuperAlpide(1), + SegmentationSuperAlpide(1), + SegmentationSuperAlpide(2), + SegmentationSuperAlpide(2)}; // TODO: fix NLayers math_utils::Point3D locCl; - if (detID >= SegmentationSuperAlpide::NLayers) { - + if (detID >= 6) { // TODO: fix NLayers o2::itsmft::SegmentationAlpide::detectorToLocalUnchecked(cl.getRow(), cl.getCol(), locCl); } else { segmentations[detID].detectorToLocalUnchecked(cl.getRow(), cl.getCol(), locCl); @@ -35,7 +40,12 @@ math_utils::Point3D TopologyDictionary::getClusterCoordinates(int detID, math_utils::Point3D TopologyDictionary::getClusterCoordinates(int detID, const its3::CompCluster& cl, const itsmft::ClusterPattern& patt, bool isGroup) { - static SegmentationSuperAlpide segmentations[SegmentationSuperAlpide::NLayers]{SegmentationSuperAlpide(0), SegmentationSuperAlpide(1), SegmentationSuperAlpide(2), SegmentationSuperAlpide(3)}; + static SegmentationSuperAlpide segmentations[6]{SegmentationSuperAlpide(0), + SegmentationSuperAlpide(0), + SegmentationSuperAlpide(1), + SegmentationSuperAlpide(1), + SegmentationSuperAlpide(2), + SegmentationSuperAlpide(2)}; // TODO: fix NLayers auto refRow = cl.getRow(); auto refCol = cl.getCol(); @@ -46,7 +56,7 @@ math_utils::Point3D TopologyDictionary::getClusterCoordinates(int detID, refCol -= round(zCOG); } math_utils::Point3D locCl; - if (detID >= SegmentationSuperAlpide::NLayers) { + if (detID >= 6) { // TODO: fix NLayers o2::itsmft::SegmentationAlpide::detectorToLocalUnchecked(refRow + xCOG, refCol + zCOG, locCl); } else { segmentations[detID].detectorToLocalUnchecked(refRow + xCOG, refCol + zCOG, locCl); diff --git a/Detectors/Upgrades/IT3/simulation/CMakeLists.txt b/Detectors/Upgrades/ITS3/simulation/CMakeLists.txt similarity index 82% rename from Detectors/Upgrades/IT3/simulation/CMakeLists.txt rename to Detectors/Upgrades/ITS3/simulation/CMakeLists.txt index 937917d4ab171..a75426ac64227 100644 --- a/Detectors/Upgrades/IT3/simulation/CMakeLists.txt +++ b/Detectors/Upgrades/ITS3/simulation/CMakeLists.txt @@ -10,14 +10,16 @@ # or submit itself to any jurisdiction. o2_add_library(ITS3Simulation - SOURCES src/ITS3Layer.cxx src/DescriptorInnerBarrelITS3.cxx + SOURCES src/ITS3Layer.cxx + src/DescriptorInnerBarrelITS3.cxx + src/Digitizer.cxx PUBLIC_LINK_LIBRARIES O2::SimulationDataFormat O2::ITSBase O2::ITSMFTSimulation ROOT::Physics) o2_target_root_dictionary(ITS3Simulation HEADERS include/ITS3Simulation/ITS3Layer.h include/ITS3Simulation/DescriptorInnerBarrelITS3.h -# include/ITS3Simulation/Digitizer.h + include/ITS3Simulation/Digitizer.h ) o2_data_file(COPY data DESTINATION Detectors/ITS3/simulation) \ No newline at end of file diff --git a/Detectors/Upgrades/IT3/simulation/data/simcuts.dat b/Detectors/Upgrades/ITS3/simulation/data/simcuts.dat similarity index 100% rename from Detectors/Upgrades/IT3/simulation/data/simcuts.dat rename to Detectors/Upgrades/ITS3/simulation/data/simcuts.dat diff --git a/Detectors/Upgrades/IT3/simulation/include/ITS3Simulation/DescriptorInnerBarrelITS3.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/DescriptorInnerBarrelITS3.h similarity index 100% rename from Detectors/Upgrades/IT3/simulation/include/ITS3Simulation/DescriptorInnerBarrelITS3.h rename to Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/DescriptorInnerBarrelITS3.h diff --git a/Detectors/Upgrades/IT3/simulation/include/ITS3Simulation/Digitizer.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h similarity index 95% rename from Detectors/Upgrades/IT3/simulation/include/ITS3Simulation/Digitizer.h rename to Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h index 1ec5fce5beb39..c933204936f8d 100644 --- a/Detectors/Upgrades/IT3/simulation/include/ITS3Simulation/Digitizer.h +++ b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h @@ -25,7 +25,7 @@ #include "ITSMFTSimulation/AlpideSimResponse.h" #include "ITSMFTSimulation/DigiParams.h" #include "ITSMFTSimulation/Hit.h" -#include "ITS3Base/GeometryTGeo.h" +#include "ITSBase/GeometryTGeo.h" #include "ITS3Base/SegmentationSuperAlpide.h" #include "DataFormatsITSMFT/Digit.h" #include "DataFormatsITSMFT/ROFRecord.h" @@ -79,7 +79,7 @@ class Digitizer : public TObject const o2::itsmft::DigiParams& getDigitParams() const { return mParams; } // provide the common itsmft::GeometryTGeo to access matrices and segmentation - void setGeometry(const o2::its3::GeometryTGeo* gm) { mGeometry = gm; } + void setGeometry(const o2::its::GeometryTGeo* gm) { mGeometry = gm; } uint32_t getEventROFrameMin() const { return mEventROFrameMin; } uint32_t getEventROFrameMax() const { return mEventROFrameMax; } @@ -107,6 +107,7 @@ class Digitizer : public TObject } std::vector mSuperSegmentations; + std::vector mLayerID; static constexpr float sec2ns = 1e9; o2::itsmft::DigiParams mParams; ///< digitization parameters @@ -122,7 +123,7 @@ class Digitizer : public TObject std::unique_ptr mAlpSimResp; // simulated response - const o2::its3::GeometryTGeo* mGeometry = nullptr; ///< ITS OR MFT upgrade geometry + const o2::its::GeometryTGeo* mGeometry = nullptr; ///< ITS OR MFT upgrade geometry std::vector mChips; ///< Array of chips digits containers std::deque> mExtraBuff; ///< burrer (per roFrame) for extra digits diff --git a/Detectors/Upgrades/IT3/simulation/include/ITS3Simulation/ITS3Layer.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ITS3Layer.h similarity index 100% rename from Detectors/Upgrades/IT3/simulation/include/ITS3Simulation/ITS3Layer.h rename to Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ITS3Layer.h diff --git a/Detectors/Upgrades/IT3/simulation/src/DescriptorInnerBarrelITS3.cxx b/Detectors/Upgrades/ITS3/simulation/src/DescriptorInnerBarrelITS3.cxx similarity index 100% rename from Detectors/Upgrades/IT3/simulation/src/DescriptorInnerBarrelITS3.cxx rename to Detectors/Upgrades/ITS3/simulation/src/DescriptorInnerBarrelITS3.cxx diff --git a/Detectors/Upgrades/IT3/simulation/src/Detector.cxx b/Detectors/Upgrades/ITS3/simulation/src/Detector.cxx similarity index 100% rename from Detectors/Upgrades/IT3/simulation/src/Detector.cxx rename to Detectors/Upgrades/ITS3/simulation/src/Detector.cxx diff --git a/Detectors/Upgrades/IT3/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx similarity index 87% rename from Detectors/Upgrades/IT3/simulation/src/Digitizer.cxx rename to Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx index 96034e8b92391..035b9b71903c2 100644 --- a/Detectors/Upgrades/IT3/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx @@ -18,6 +18,7 @@ #include "MathUtils/Cartesian.h" #include "SimulationDataFormat/MCTruthContainer.h" #include "DetectorsRaw/HBFUtils.h" +#include "CommonConstants/MathConstants.h" #include #include @@ -37,11 +38,17 @@ using namespace o2::its3; void Digitizer::init() { - for (int iLayer{0}; iLayer < SegmentationSuperAlpide::NLayers; ++iLayer) { - mSuperSegmentations.push_back(SegmentationSuperAlpide(iLayer)); + mLayerID.clear(); + mSuperSegmentations.clear(); + for (int iLayer{0}; iLayer < mGeometry->getNumberOfLayers() - 4; ++iLayer) { + for (int iChip{0}; iChip < mGeometry->getNumberOfChipsPerLayer(iLayer); ++iChip) { + LOGP(info, "layer: {} chip: {}", iLayer, iChip); + mSuperSegmentations.push_back(SegmentationSuperAlpide(iLayer)); + mLayerID.push_back(iLayer); + } } - const Int_t numOfChips = mGeometry->getNumberOfChips() + SegmentationSuperAlpide::NLayers; + const Int_t numOfChips = mGeometry->getNumberOfChips(); mChips.resize(numOfChips); for (int i = numOfChips; i--;) { mChips[i].setChipIndex(i); @@ -142,7 +149,7 @@ void Digitizer::fillOutputContainer(uint32_t frameLast) auto& extra = *(mExtraBuff.front().get()); for (int iChip{0}; iChip < mChips.size(); ++iChip) { auto& chip = mChips[iChip]; - if (iChip < SegmentationSuperAlpide::NLayers) { + if (iChip < mSuperSegmentations.size()) { // Check if this is a chip of ITS3 chip.addNoise(mROFrameMin, mROFrameMin, &mParams, mSuperSegmentations[iChip].NRows, mSuperSegmentations[iChip].NCols); } else { chip.addNoise(mROFrameMin, mROFrameMin, &mParams); @@ -224,16 +231,28 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID float nStepsInv = mParams.getNSimStepsInv(); int nSteps = mParams.getNSimSteps(); short detID{hit.GetDetectorID()}; - const auto& matrix = mGeometry->getMatrixL2G(detID); - bool innerBarrel{detID < SegmentationSuperAlpide::NLayers}; + const auto& matrix = mGeometry->getMatrixL2G(detID); // <<<< ????? + bool innerBarrel{detID < mSuperSegmentations.size()}; + float gap = 0.1; // FIXME: get this properly! auto startPos = hit.GetPosStart(); - float startPhi{std::atan2(-startPos.Y(), -startPos.X())}; - auto endPos = hit.GetPos(); - float endPhi{std::atan2(-endPos.Y(), -endPos.X())}; + // LOGP(info, "StartPos: {} {} {}", startPos.X(), startPos.Y(), startPos.Z()); math_utils::Vector3D xyzLocS, xyzLocE; if (innerBarrel) { - xyzLocS = {SegmentationSuperAlpide::Radii[detID] * startPhi, 0.f, startPos.Z()}; - xyzLocE = {SegmentationSuperAlpide::Radii[detID] * endPhi, 0.f, endPos.Z()}; + bool isTop = startPos.Y() > 0; + float reShiftedStartY = isTop ? startPos.Y() - gap / 2 : startPos.Y() + gap / 2; // This is due to the gap between the ITS3 emispheres + float startPhi{std::atan2(-reShiftedStartY, -startPos.X()) + (isTop ? constants::math::PI / 2 : -constants::math::PI / 2)}; + if (innerBarrel) { + // LOGP(info, "X: {}, Y: {}, startPhi: {}", startPos.X(), reShiftedStartY, startPhi); + } + auto endPos = hit.GetPos(); + float reShiftedEndY = isTop ? endPos.Y() - gap / 2 : endPos.Y() + gap / 2; // This is due to the gap between the ITS3 emispheres + float endPhi{std::atan2(-reShiftedEndY, -endPos.X()) + (isTop ? constants::math::PI / 2 : -constants::math::PI / 2)}; + if (innerBarrel) { + // LOGP(info, "X: {}, Y: {}, endPhi: {}", endPos.X(), reShiftedEndY, endPhi); + } + float deltaY = reShiftedEndY - reShiftedEndY; + xyzLocS = {SegmentationSuperAlpide::Radii[mLayerID[detID]] * (startPhi), 0.f, startPos.Z()}; + xyzLocE = {SegmentationSuperAlpide::Radii[mLayerID[detID]] * (endPhi), deltaY, endPos.Z()}; } else { xyzLocS = matrix ^ (hit.GetPosStart()); xyzLocE = matrix ^ (hit.GetPos()); @@ -242,16 +261,17 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID math_utils::Vector3D step(xyzLocE); step -= xyzLocS; step *= nStepsInv; // position increment at each step - // the electrons will injected in the middle of each step + // the electrons will be injected in the middle of each step math_utils::Vector3D stepH(step * 0.5); - xyzLocS += stepH; - xyzLocE -= stepH; + xyzLocS += stepH; // Adjust start position to the middle of the first step + xyzLocE -= stepH; // Adjust end position to the middle of the last step int rowS = -1, colS = -1, rowE = -1, colE = -1, nSkip = 0; if (innerBarrel) { // get entrance pixel row and col while (!mSuperSegmentations[detID].localToDetector(xyzLocS.X(), xyzLocS.Z(), rowS, colS)) { // guard-ring ? if (++nSkip >= nSteps) { + LOGP(info, "Start: detId {}", detID); return; // did not enter to sensitive matrix } xyzLocS += step; @@ -259,10 +279,12 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID // get exit pixel row and col while (!mSuperSegmentations[detID].localToDetector(xyzLocE.X(), xyzLocE.Z(), rowE, colE)) { // guard-ring ? if (++nSkip >= nSteps) { + LOGP(info, "End: detId {}", detID); return; // did not enter to sensitive matrix } xyzLocE -= step; } + // LOGP(info, "x: {}, z: {}, col: {}, row: {}, detID: {}", xyzLocS.X(), xyzLocS.Z(), colS, rowS, detID); } else { // get entrance pixel row and col while (!Segmentation::localToDetector(xyzLocS.X(), xyzLocS.Z(), rowS, colS)) { // guard-ring ? diff --git a/Detectors/Upgrades/IT3/simulation/src/ITS3Layer.cxx b/Detectors/Upgrades/ITS3/simulation/src/ITS3Layer.cxx similarity index 100% rename from Detectors/Upgrades/IT3/simulation/src/ITS3Layer.cxx rename to Detectors/Upgrades/ITS3/simulation/src/ITS3Layer.cxx diff --git a/Detectors/Upgrades/IT3/simulation/src/ITS3SimulationLinkDef.h b/Detectors/Upgrades/ITS3/simulation/src/ITS3SimulationLinkDef.h similarity index 94% rename from Detectors/Upgrades/IT3/simulation/src/ITS3SimulationLinkDef.h rename to Detectors/Upgrades/ITS3/simulation/src/ITS3SimulationLinkDef.h index 7a3dff2e9026b..bb9c93b6e20e6 100644 --- a/Detectors/Upgrades/IT3/simulation/src/ITS3SimulationLinkDef.h +++ b/Detectors/Upgrades/ITS3/simulation/src/ITS3SimulationLinkDef.h @@ -17,6 +17,6 @@ #pragma link C++ class o2::its3::ITS3Layer + ; #pragma link C++ class o2::its3::DescriptorInnerBarrelITS3 + ; -//#pragma link C++ class o2::its3::Digitizer + ; +#pragma link C++ class o2::its3::Digitizer + ; #endif diff --git a/Detectors/Upgrades/IT3/workflow/CMakeLists.txt b/Detectors/Upgrades/ITS3/workflow/CMakeLists.txt similarity index 79% rename from Detectors/Upgrades/IT3/workflow/CMakeLists.txt rename to Detectors/Upgrades/ITS3/workflow/CMakeLists.txt index 7aca1a47c7743..ec3ec7ab390e6 100644 --- a/Detectors/Upgrades/IT3/workflow/CMakeLists.txt +++ b/Detectors/Upgrades/ITS3/workflow/CMakeLists.txt @@ -13,25 +13,26 @@ o2_add_library(ITS3Workflow SOURCES src/DigitReaderSpec.cxx src/DigitWriterSpec.cxx src/RecoWorkflow.cxx - src/ClusterWriterWorkflow.cxx + # src/ClusterWriterWorkflow.cxx src/ClustererSpec.cxx src/ClusterWriterSpec.cxx - # src/TrackerSpec.cxx + src/TrackerSpec.cxx # src/CookedTrackerSpec.cxx # src/TrackWriterSpec.cxx # src/TrackReaderSpec.cxx # src/VertexReaderSpec.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::SimConfig - O2::DataFormatsITS + O2::DataFormatsITSMFT O2::DataFormatsITS3 O2::SimulationDataFormat O2::ITS3Simulation - # O2::ITStracking - # O2::ITSReconstruction + O2::ITStracking + O2::ITSMFTReconstruction O2::ITS3Reconstruction - # O2::ITSMFTWorkflow + O2::ITSWorkflow O2::GPUTracking + O2::ITSBase ) # o2_add_executable(digit-writer-workflow @@ -44,10 +45,10 @@ o2_add_library(ITS3Workflow # COMPONENT_NAME its3 # PUBLIC_LINK_LIBRARIES O2::ITS3Workflow) -# o2_add_executable(reco-workflow -# SOURCES src/its3-reco-workflow.cxx -# COMPONENT_NAME its3 -# PUBLIC_LINK_LIBRARIES O2::ITS3Workflow) +o2_add_executable(reco-workflow + SOURCES src/its3-reco-workflow.cxx + COMPONENT_NAME its3 + PUBLIC_LINK_LIBRARIES O2::ITS3Workflow) # o2_add_executable(cluster-writer-workflow # SOURCES src/its-cluster-writer-workflow.cxx diff --git a/Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/ClusterReaderSpec.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClusterReaderSpec.h similarity index 100% rename from Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/ClusterReaderSpec.h rename to Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClusterReaderSpec.h diff --git a/Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/ClusterWriterSpec.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClusterWriterSpec.h similarity index 100% rename from Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/ClusterWriterSpec.h rename to Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClusterWriterSpec.h diff --git a/Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/ClusterWriterWorkflow.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClusterWriterWorkflow.h similarity index 100% rename from Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/ClusterWriterWorkflow.h rename to Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClusterWriterWorkflow.h diff --git a/Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/ClustererSpec.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h similarity index 95% rename from Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/ClustererSpec.h rename to Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h index c1dd586a86766..3a96921812c67 100644 --- a/Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/ClustererSpec.h +++ b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h @@ -11,8 +11,8 @@ /// @file ClustererSpec.h -#ifndef O2_ITS_CLUSTERERDPL -#define O2_ITS_CLUSTERERDPL +#ifndef O2_ITS3_CLUSTERERDPL +#define O2_ITS3_CLUSTERERDPL #include diff --git a/Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/DigitReaderSpec.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitReaderSpec.h similarity index 100% rename from Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/DigitReaderSpec.h rename to Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitReaderSpec.h diff --git a/Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/DigitWriterSpec.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitWriterSpec.h similarity index 100% rename from Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/DigitWriterSpec.h rename to Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitWriterSpec.h diff --git a/Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/RecoWorkflow.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/RecoWorkflow.h similarity index 62% rename from Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/RecoWorkflow.h rename to Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/RecoWorkflow.h index 0d4dbe0564e41..eeeb2d3742df2 100644 --- a/Detectors/Upgrades/IT3/workflow/include/ITS3Workflow/RecoWorkflow.h +++ b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/RecoWorkflow.h @@ -9,8 +9,8 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -#ifndef O2_ITS_RECOWORKFLOW_H -#define O2_ITS_RECOWORKFLOW_H +#ifndef O2_ITS3_RECOWORKFLOW_H +#define O2_ITS3_RECOWORKFLOW_H /// @file RecoWorkflow.h @@ -28,8 +28,13 @@ namespace its3 namespace reco_workflow { -framework::WorkflowSpec getWorkflow(bool useMC, const std::string& trmode, o2::gpu::GPUDataTypes::DeviceType dType = o2::gpu::GPUDataTypes::DeviceType::CPU, - bool upstreamDigits = false, bool upstreamClusters = false, bool disableRootOutput = false, bool eencode = false); +framework::WorkflowSpec getWorkflow(bool useMC, + const std::string& trmode, + o2::gpu::GPUDataTypes::DeviceType dtype, + bool upstreamDigits, + bool upstreamClusters, + bool disableRootOutput, + int useTrig); } } // namespace its3 diff --git a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/TrackerSpec.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/TrackerSpec.h new file mode 100644 index 0000000000000..64e9f86ad4142 --- /dev/null +++ b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/TrackerSpec.h @@ -0,0 +1,76 @@ +// 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 TrackerSpec.h + +#ifndef O2_ITS_TRACKERDPL +#define O2_ITS_TRACKERDPL + +#include "DataFormatsParameters/GRPObject.h" +#include "DataFormatsITSMFT/TopologyDictionary.h" + +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" + +#include "ITStracking/TimeFrame.h" +#include "ITStracking/Tracker.h" +#include "ITStracking/TrackerTraits.h" +#include "ITStracking/Vertexer.h" +#include "ITStracking/VertexerTraits.h" + +#include "GPUO2Interface.h" +#include "GPUReconstruction.h" +#include "GPUChainITS.h" +#include "CommonUtils/StringUtils.h" +#include "TStopwatch.h" +#include "DetectorsBase/GRPGeomHelper.h" + +namespace o2 +{ +namespace its3 +{ + +class TrackerDPL : public framework::Task +{ + public: + TrackerDPL(std::shared_ptr gr, bool isMC, int trgType, const std::string& trModeS, o2::gpu::GPUDataTypes::DeviceType dType = o2::gpu::GPUDataTypes::DeviceType::CPU); + ~TrackerDPL() override = default; + void init(framework::InitContext& ic) final; + void run(framework::ProcessingContext& pc) final; + void endOfStream(framework::EndOfStreamContext& ec) final; + void finaliseCCDB(framework::ConcreteDataMatcher& matcher, void* obj) final; + void setClusterDictionary(const o2::itsmft::TopologyDictionary* d) { mDict = d; } + + private: + void updateTimeDependentParams(framework::ProcessingContext& pc); + + bool mIsMC = false; + bool mRunVertexer = true; + bool mCosmicsProcessing = false; + int mUseTriggers = 0; + std::string mMode = "sync"; + std::shared_ptr mGGCCDBRequest; + const o2::itsmft::TopologyDictionary* mDict = nullptr; + std::unique_ptr mRecChain = nullptr; + std::unique_ptr mChainITS = nullptr; + std::unique_ptr mTracker = nullptr; + std::unique_ptr mVertexer = nullptr; + TStopwatch mTimer; +}; + +/// create a processor spec +/// run ITS CA tracker +framework::DataProcessorSpec getTrackerSpec(bool useMC, int useTrig, const std::string& trModeS, o2::gpu::GPUDataTypes::DeviceType dType); + +} // namespace its3 +} // namespace o2 + +#endif /* O2_ITS_TRACKERDPL */ diff --git a/Detectors/Upgrades/IT3/workflow/src/ClusterReaderSpec.cxx b/Detectors/Upgrades/ITS3/workflow/src/ClusterReaderSpec.cxx similarity index 100% rename from Detectors/Upgrades/IT3/workflow/src/ClusterReaderSpec.cxx rename to Detectors/Upgrades/ITS3/workflow/src/ClusterReaderSpec.cxx diff --git a/Detectors/Upgrades/IT3/workflow/src/ClusterWriterSpec.cxx b/Detectors/Upgrades/ITS3/workflow/src/ClusterWriterSpec.cxx similarity index 100% rename from Detectors/Upgrades/IT3/workflow/src/ClusterWriterSpec.cxx rename to Detectors/Upgrades/ITS3/workflow/src/ClusterWriterSpec.cxx diff --git a/Detectors/Upgrades/IT3/workflow/src/ClusterWriterWorkflow.cxx b/Detectors/Upgrades/ITS3/workflow/src/ClusterWriterWorkflow.cxx similarity index 100% rename from Detectors/Upgrades/IT3/workflow/src/ClusterWriterWorkflow.cxx rename to Detectors/Upgrades/ITS3/workflow/src/ClusterWriterWorkflow.cxx diff --git a/Detectors/Upgrades/IT3/workflow/src/ClustererSpec.cxx b/Detectors/Upgrades/ITS3/workflow/src/ClustererSpec.cxx similarity index 98% rename from Detectors/Upgrades/IT3/workflow/src/ClustererSpec.cxx rename to Detectors/Upgrades/ITS3/workflow/src/ClustererSpec.cxx index 93890fa01fffe..44951be8a442d 100644 --- a/Detectors/Upgrades/IT3/workflow/src/ClustererSpec.cxx +++ b/Detectors/Upgrades/ITS3/workflow/src/ClustererSpec.cxx @@ -40,7 +40,7 @@ namespace its3 void ClustererDPL::init(InitContext& ic) { mClusterer = std::make_unique(); - mClusterer->setNChips(o2::itsmft::ChipMappingITS::getNChips(o2::itsmft::ChipMappingITS::MB) + o2::itsmft::ChipMappingITS::getNChips(o2::itsmft::ChipMappingITS::OB) + SegmentationSuperAlpide::NLayers); + mClusterer->setNChips(o2::itsmft::ChipMappingITS::getNChips(o2::itsmft::ChipMappingITS::MB) + o2::itsmft::ChipMappingITS::getNChips(o2::itsmft::ChipMappingITS::OB) + 6); /// Fix MEEEE auto filenameGRP = ic.options().get("grp-file"); const auto grp = o2::parameters::GRPObject::loadFrom(filenameGRP.c_str()); diff --git a/Detectors/Upgrades/IT3/workflow/src/DigitReaderSpec.cxx b/Detectors/Upgrades/ITS3/workflow/src/DigitReaderSpec.cxx similarity index 100% rename from Detectors/Upgrades/IT3/workflow/src/DigitReaderSpec.cxx rename to Detectors/Upgrades/ITS3/workflow/src/DigitReaderSpec.cxx diff --git a/Detectors/Upgrades/IT3/workflow/src/DigitWriterSpec.cxx b/Detectors/Upgrades/ITS3/workflow/src/DigitWriterSpec.cxx similarity index 100% rename from Detectors/Upgrades/IT3/workflow/src/DigitWriterSpec.cxx rename to Detectors/Upgrades/ITS3/workflow/src/DigitWriterSpec.cxx diff --git a/Detectors/Upgrades/IT3/workflow/src/RecoWorkflow.cxx b/Detectors/Upgrades/ITS3/workflow/src/RecoWorkflow.cxx similarity index 73% rename from Detectors/Upgrades/IT3/workflow/src/RecoWorkflow.cxx rename to Detectors/Upgrades/ITS3/workflow/src/RecoWorkflow.cxx index 325687aafac48..653752271cda7 100644 --- a/Detectors/Upgrades/IT3/workflow/src/RecoWorkflow.cxx +++ b/Detectors/Upgrades/ITS3/workflow/src/RecoWorkflow.cxx @@ -9,14 +9,11 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// @file RecoWorkflow.cxx - #include "ITS3Workflow/RecoWorkflow.h" #include "ITS3Workflow/ClustererSpec.h" #include "ITS3Workflow/ClusterWriterSpec.h" -// #include "ITSWorkflow/TrackerSpec.h" -// #include "ITSWorkflow/TrackWriterSpec.h" -// #include "ITSMFTWorkflow/EntropyEncoderSpec.h" +#include "ITSWorkflow/TrackerSpec.h" +#include "ITSWorkflow/TrackWriterSpec.h" #include "ITS3Workflow/DigitReaderSpec.h" namespace o2 @@ -28,8 +25,7 @@ namespace reco_workflow { framework::WorkflowSpec getWorkflow(bool useMC, const std::string& trmode, o2::gpu::GPUDataTypes::DeviceType dtype, - bool upstreamDigits, bool upstreamClusters, bool disableRootOutput, - bool) + bool upstreamDigits, bool upstreamClusters, bool disableRootOutput, int useTrig) { framework::WorkflowSpec specs; @@ -40,17 +36,16 @@ framework::WorkflowSpec getWorkflow(bool useMC, const std::string& trmode, o2::g if (!upstreamClusters) { specs.emplace_back(o2::its3::getClustererSpec(useMC)); } + if (!disableRootOutput) { specs.emplace_back(o2::its3::getClusterWriterSpec(useMC)); } - // specs.emplace_back(o2::its::getTrackerSpec(useMC, trmode, dtype)); - // if (!disableRootOutput) { - // specs.emplace_back(o2::its::getTrackWriterSpec(useMC)); - // } - - // if (eencode) { - // specs.emplace_back(o2::itsmft::getEntropyEncoderSpec("ITS")); - // } + + specs.emplace_back(o2::its::getTrackerSpec(useMC, useTrig, trmode, dtype)); + if (!disableRootOutput) { + specs.emplace_back(o2::its::getTrackWriterSpec(useMC)); + } + return specs; } diff --git a/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx b/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx new file mode 100644 index 0000000000000..d854f136a75cc --- /dev/null +++ b/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx @@ -0,0 +1,414 @@ +// 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 TrackerSpec.cxx + +#include + +#include "TGeoGlobalMagField.h" + +#include "Framework/ControlService.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/CCDBParamSpec.h" +#include "ITS3Workflow/TrackerSpec.h" +#include "DataFormatsITS3/CompCluster.h" +#include "DataFormatsITS/TrackITS.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "DataFormatsITSMFT/PhysTrigger.h" + +#include "ITStracking/ROframe.h" +#include "ITStracking/IOUtils.h" +#include "ITStracking/TrackingConfigParam.h" +#include "ITSMFTBase/DPLAlpideParam.h" +#include "ITSMFTReconstruction/ClustererParam.h" + +#include "Field/MagneticField.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" +#include "ITSBase/GeometryTGeo.h" +#include "CommonDataFormat/IRFrame.h" +#include "DetectorsCommonDataFormats/DetectorNameConf.h" +#include "DataFormatsTRD/TriggerRecord.h" +#include "ITSReconstruction/FastMultEstConfig.h" +#include "ITSReconstruction/FastMultEst.h" +#include + +namespace o2 +{ +using namespace framework; +using its::FastMultEst; +using its::FastMultEstConfig; +using its::TimeFrame; +using its::Tracker; +using its::TrackingParameters; +using its::TrackITSExt; +using its::Vertexer; + +namespace its3 +{ +using Vertex = o2::dataformats::Vertex>; + +TrackerDPL::TrackerDPL(std::shared_ptr gr, bool isMC, int trgType, const std::string& trModeS, o2::gpu::GPUDataTypes::DeviceType dType) : mGGCCDBRequest(gr), mIsMC{isMC}, mUseTriggers{trgType}, mMode{trModeS}, mRecChain{o2::gpu::GPUReconstruction::CreateInstance(dType, true)} +{ + std::transform(mMode.begin(), mMode.end(), mMode.begin(), [](unsigned char c) { return std::tolower(c); }); +} + +void TrackerDPL::init(InitContext& ic) +{ + mTimer.Stop(); + mTimer.Reset(); + o2::base::GRPGeomHelper::instance().setRequest(mGGCCDBRequest); + mChainITS.reset(mRecChain->AddChain()); + mVertexer = std::make_unique(mChainITS->GetITSVertexerTraits()); + mTracker = std::make_unique(mChainITS->GetITSTrackerTraits()); + mRunVertexer = true; + mCosmicsProcessing = false; + std::vector trackParams; + + if (mMode == "async") { + + trackParams.resize(3); + trackParams[1].TrackletMinPt = 0.2f; + trackParams[1].CellDeltaTanLambdaSigma *= 2.; + trackParams[2].TrackletMinPt = 0.1f; + trackParams[2].CellDeltaTanLambdaSigma *= 4.; + trackParams[2].MinTrackLength = 4; + LOG(info) << "Initializing tracker in async. phase reconstruction with " << trackParams.size() << " passes"; + + } else if (mMode == "sync_misaligned") { + + trackParams.resize(3); + trackParams[0].PhiBins = 32; + trackParams[0].ZBins = 64; + trackParams[0].CellDeltaTanLambdaSigma *= 3.; + trackParams[0].LayerMisalignment[0] = 1.e-2; + trackParams[0].LayerMisalignment[1] = 1.e-2; + trackParams[0].LayerMisalignment[2] = 1.e-2; + trackParams[0].LayerMisalignment[3] = 3.e-2; + trackParams[0].LayerMisalignment[4] = 3.e-2; + trackParams[0].LayerMisalignment[5] = 3.e-2; + trackParams[0].LayerMisalignment[6] = 3.e-2; + trackParams[0].FitIterationMaxChi2[0] = 50.; + trackParams[0].FitIterationMaxChi2[1] = 25.; + trackParams[1] = trackParams[0]; + trackParams[2] = trackParams[0]; + trackParams[1].MinTrackLength = 6; + trackParams[2].MinTrackLength = 4; + LOG(info) << "Initializing tracker in misaligned sync. phase reconstruction with " << trackParams.size() << " passes"; + + } else if (mMode == "sync") { + trackParams.resize(1); + LOG(info) << "Initializing tracker in sync. phase reconstruction with " << trackParams.size() << " passes"; + } else if (mMode == "cosmics") { + mCosmicsProcessing = true; + mRunVertexer = false; + trackParams.resize(1); + trackParams[0].MinTrackLength = 4; + trackParams[0].CellDeltaTanLambdaSigma *= 10; + trackParams[0].PhiBins = 4; + trackParams[0].ZBins = 16; + trackParams[0].PVres = 1.e5f; + trackParams[0].LayerMisalignment[0] = 1.e-2; + trackParams[0].LayerMisalignment[1] = 1.e-2; + trackParams[0].LayerMisalignment[2] = 1.e-2; + trackParams[0].LayerMisalignment[3] = 3.e-2; + trackParams[0].LayerMisalignment[4] = 3.e-2; + trackParams[0].LayerMisalignment[5] = 3.e-2; + trackParams[0].LayerMisalignment[6] = 3.e-2; + trackParams[0].FitIterationMaxChi2[0] = 50.; + trackParams[0].FitIterationMaxChi2[1] = 25.; + trackParams[0].TrackletsPerClusterLimit = 100.; + trackParams[0].CellsPerClusterLimit = 100.; + LOG(info) << "Initializing tracker in reconstruction for cosmics with " << trackParams.size() << " passes"; + + } else { + throw std::runtime_error(fmt::format("Unsupported ITS3 tracking mode {:s} ", mMode)); + } + + for (auto& params : trackParams) { + params.CorrType = o2::base::PropagatorImpl::MatCorrType::USEMatCorrLUT; + } + mTracker->setParameters(trackParams); +} + +void TrackerDPL::run(ProcessingContext& pc) +{ + mTimer.Start(false); + updateTimeDependentParams(pc); + auto compClusters = pc.inputs().get>("compClusters"); + gsl::span patterns = pc.inputs().get>("patterns"); + gsl::span physTriggers; + std::vector fromTRD; + if (mUseTriggers == 2) { // use TRD triggers + o2::InteractionRecord ir{0, pc.services().get().firstTForbit}; + auto trdTriggers = pc.inputs().get>("phystrig"); + for (const auto& trig : trdTriggers) { + if (trig.getBCData() >= ir && trig.getNumberOfTracklets()) { + ir = trig.getBCData(); + fromTRD.emplace_back(o2::itsmft::PhysTrigger{ir, 0}); + } + } + physTriggers = gsl::span(fromTRD.data(), fromTRD.size()); + } else if (mUseTriggers == 1) { // use Phys triggers from ITS stream + physTriggers = pc.inputs().get>("phystrig"); + } + + // code further down does assignment to the rofs and the altered object is used for output + // we therefore need a copy of the vector rather than an object created directly on the input data, + // the output vector however is created directly inside the message memory thus avoiding copy by + // snapshot + auto rofsinput = pc.inputs().get>("ROframes"); + auto& rofs = pc.outputs().make>(Output{"ITS3", "ITS3TrackROF", 0, Lifetime::Timeframe}, rofsinput.begin(), rofsinput.end()); + + auto& irFrames = pc.outputs().make>(Output{"ITS3", "IRFRAMES", 0, Lifetime::Timeframe}); + + const auto& alpParams = o2::itsmft::DPLAlpideParam::Instance(); // RS: this should come from CCDB + int nBCPerTF = alpParams.roFrameLengthInBC; + + LOG(info) << "ITS3Tracker pulled " << compClusters.size() << " clusters, " << rofs.size() << " RO frames"; + + const dataformats::MCTruthContainer* labels = nullptr; + gsl::span mc2rofs; + if (mIsMC) { + labels = pc.inputs().get*>("labels").release(); + // get the array as read-only span, a snapshot is send forward + mc2rofs = pc.inputs().get>("MC2ROframes"); + LOG(info) << labels->getIndexedSize() << " MC label objects , in " << mc2rofs.size() << " MC events"; + } + + std::vector tracks; + auto& allClusIdx = pc.outputs().make>(Output{"ITS3", "TRACKCLSID", 0, Lifetime::Timeframe}); + std::vector trackLabels; + std::vector verticesLabels; + auto& allTracks = pc.outputs().make>(Output{"ITS3", "TRACKS", 0, Lifetime::Timeframe}); + std::vector allTrackLabels; + std::vector allVerticesLabels; + + auto& vertROFvec = pc.outputs().make>(Output{"ITS3", "VERTICESROF", 0, Lifetime::Timeframe}); + auto& vertices = pc.outputs().make>(Output{"ITS3", "VERTICES", 0, Lifetime::Timeframe}); + + std::uint32_t roFrame = 0; + + bool continuous = o2::base::GRPGeomHelper::instance().getGRPECS()->isDetContinuousReadOut(o2::detectors::DetID::IT3); + LOG(info) << "ITS3Tracker RO: continuous=" << continuous; + TimeFrame* timeFrame = mChainITS->GetITSTimeframe(); + mTracker->adoptTimeFrame(*timeFrame); + + mTracker->setBz(o2::base::Propagator::Instance()->getNominalBz()); + mVertexer->adoptTimeFrame(*timeFrame); + gsl::span::iterator pattIt = patterns.begin(); + + gsl::span rofspan(rofs); + timeFrame->loadROFrameDataITS3(rofspan, compClusters, pattIt, mDict, labels); + pattIt = patterns.begin(); + std::vector savedROF; + auto logger = [&](std::string s) { LOG(info) << s; }; + auto errorLogger = [&](std::string s) { LOG(error) << s; }; + + FastMultEst multEst; // mult estimator + std::vector processingMask; + int cutVertexMult{0}, cutRandomMult = int(rofs.size()) - multEst.selectROFs(rofs, compClusters, physTriggers, processingMask); + timeFrame->setMultiplicityCutMask(processingMask); + float vertexerElapsedTime{0.f}; + if (mRunVertexer) { + // Run seeding vertexer + vertexerElapsedTime = mVertexer->clustersToVertices(logger); + } + const auto& multEstConf = FastMultEstConfig::Instance(); // parameters for mult estimation and cuts + for (auto iRof{0}; iRof < rofspan.size(); ++iRof) { + std::vector vtxVecLoc; + auto& vtxROF = vertROFvec.emplace_back(rofspan[iRof]); + vtxROF.setFirstEntry(vertices.size()); + if (mRunVertexer) { + auto vtxSpan = timeFrame->getPrimaryVertices(iRof); + vtxROF.setNEntries(vtxSpan.size()); + bool selROF = vtxSpan.size() == 0; + for (auto iV{0}; iV < vtxSpan.size(); ++iV) { + auto& v = vtxSpan[iV]; + if (multEstConf.isVtxMultCutRequested() && !multEstConf.isPassingVtxMultCut(v.getNContributors())) { + continue; // skip vertex of unwanted multiplicity + } + selROF = true; + vertices.push_back(v); + if (mIsMC) { + auto vLabels = timeFrame->getPrimaryVerticesLabels(iRof)[iV]; + std::copy(vLabels.begin(), vLabels.end(), std::back_inserter(allVerticesLabels)); + } + } + if (processingMask[iRof] && !selROF) { // passed selection in clusters and not in vertex multiplicity + LOG(debug) << fmt::format("ROF {} rejected by the vertex multiplicity selection [{},{}]", + iRof, + multEstConf.cutMultVtxLow, + multEstConf.cutMultVtxHigh); + processingMask[iRof] = selROF; + cutVertexMult++; + } + } else { // cosmics + vtxVecLoc.emplace_back(Vertex()); + vtxVecLoc.back().setNContributors(1); + vtxROF.setNEntries(vtxVecLoc.size()); + for (auto& v : vtxVecLoc) { + vertices.push_back(v); + } + timeFrame->addPrimaryVertices(vtxVecLoc); + } + } + LOG(info) << fmt::format(" - rejected {}/{} ROFs: random/mult.sel:{} (seed {}), vtx.sel:{}", cutRandomMult + cutVertexMult, rofspan.size(), cutRandomMult, multEst.lastRandomSeed, cutVertexMult); + LOG(info) << fmt::format(" - Vertex seeding total elapsed time: {} ms for {} vertices found in {} ROFs", vertexerElapsedTime, timeFrame->getPrimaryVerticesNum(), rofspan.size()); + LOG(info) << fmt::format(" - Beam position computed for the TF: {}, {}", timeFrame->getBeamX(), timeFrame->getBeamY()); + + if (mCosmicsProcessing && compClusters.size() > 1500 * rofspan.size()) { + LOG(error) << "Cosmics processing was requested with an average detector occupancy exceeding 1.e-7, skipping TF processing."; + } else { + + timeFrame->setMultiplicityCutMask(processingMask); + // Run CA tracker + mTracker->clustersToTracks(logger, errorLogger); + if (timeFrame->hasBogusClusters()) { + LOG(warning) << fmt::format(" - The processed timeframe had {} clusters with wild z coordinates, check the dictionaries", timeFrame->hasBogusClusters()); + } + + for (unsigned int iROF{0}; iROF < rofs.size(); ++iROF) { + auto& rof{rofs[iROF]}; + tracks = timeFrame->getTracks(iROF); + trackLabels = timeFrame->getTracksLabel(iROF); + auto number{tracks.size()}; + auto first{allTracks.size()}; + int offset = -rof.getFirstEntry(); // cluster entry!!! + rof.setFirstEntry(first); + rof.setNEntries(number); + + if (processingMask[iROF]) { + irFrames.emplace_back(rof.getBCData(), rof.getBCData() + nBCPerTF - 1).info = tracks.size(); + } + + 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; + 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++; + } + } + assert(ncl == nclf); + allTracks.emplace_back(trc); + } + } + LOGP(info, "ITS3Tracker pushed {} tracks and {} vertices", allTracks.size(), vertices.size()); + if (mIsMC) { + LOGP(info, "ITS3Tracker pushed {} track labels", allTrackLabels.size()); + LOGP(info, "ITS3Tracker pushed {} vertex labels", allVerticesLabels.size()); + + pc.outputs().snapshot(Output{"ITS3", "TRACKSMCTR", 0, Lifetime::Timeframe}, allTrackLabels); + pc.outputs().snapshot(Output{"ITS3", "VERTICESMCTR", 0, Lifetime::Timeframe}, allVerticesLabels); + pc.outputs().snapshot(Output{"ITS3", "ITSTrackMC2ROF", 0, Lifetime::Timeframe}, mc2rofs); + } + } + mTimer.Stop(); +} + +///_______________________________________ +void TrackerDPL::updateTimeDependentParams(ProcessingContext& pc) +{ + o2::base::GRPGeomHelper::instance().checkUpdates(pc); + static bool initOnceDone = false; + if (!initOnceDone) { // this params need to be queried only once + initOnceDone = true; + pc.inputs().get("cldict"); // just to trigger the finaliseCCDB + pc.inputs().get*>("alppar"); + GeometryTGeo* geom = GeometryTGeo::Instance(); + geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::T2GRot, o2::math_utils::TransformType::T2G)); + mVertexer->getGlobalConfiguration(); + mTracker->getGlobalConfiguration(); + } +} + +///_______________________________________ +void TrackerDPL::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) +{ + if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { + return; + } + if (matcher == ConcreteDataMatcher("ITS", "CLUSDICT", 0)) { + LOG(info) << "cluster dictionary updated"; + setClusterDictionary((const o2::itsmft::TopologyDictionary*)obj); + return; + } + // Note: strictly speaking, for Configurable params we don't need finaliseCCDB check, the singletons are updated at the CCDB fetcher level + if (matcher == ConcreteDataMatcher("ITS", "ALPIDEPARAM", 0)) { + LOG(info) << "Alpide param updated"; + const auto& par = o2::itsmft::DPLAlpideParam::Instance(); + par.printKeyValues(); + return; + } +} + +void TrackerDPL::endOfStream(EndOfStreamContext& ec) +{ + LOGF(info, "ITS3 CA-Tracker total timing: Cpu: %.3e Real: %.3e s in %d slots", + mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); +} + +DataProcessorSpec getTrackerSpec(bool useMC, int trgType, const std::string& trModeS, o2::gpu::GPUDataTypes::DeviceType dType) +{ + std::vector inputs; + inputs.emplace_back("compClusters", "ITS3", "COMPCLUSTERS", 0, Lifetime::Timeframe); + inputs.emplace_back("patterns", "ITS3", "PATTERNS", 0, Lifetime::Timeframe); + inputs.emplace_back("ROframes", "ITS3", "CLUSTERSROF", 0, Lifetime::Timeframe); + if (trgType == 1) { + inputs.emplace_back("phystrig", "ITS3", "PHYSTRIG", 0, Lifetime::Timeframe); + } else if (trgType == 2) { + inputs.emplace_back("phystrig", "TRD", "TRKTRGRD", 0, Lifetime::Timeframe); + } + inputs.emplace_back("cldict", "ITS", "CLUSDICT", 0, Lifetime::Condition, ccdbParamSpec("ITS/Calib/ClusterDictionary")); + inputs.emplace_back("alppar", "ITS", "ALPIDEPARAM", 0, Lifetime::Condition, ccdbParamSpec("ITS/Config/AlpideParam")); + auto ggRequest = std::make_shared(false, // orbitResetTime + true, // GRPECS=true + false, // GRPLHCIF + true, // GRPMagField + true, // askMatLUT + o2::base::GRPGeomRequest::Aligned, // geometry + inputs, + true); + std::vector outputs; + outputs.emplace_back("ITS3", "TRACKS", 0, Lifetime::Timeframe); + outputs.emplace_back("ITS3", "TRACKCLSID", 0, Lifetime::Timeframe); + outputs.emplace_back("ITS3", "ITSTrackROF", 0, Lifetime::Timeframe); + outputs.emplace_back("ITS3", "VERTICES", 0, Lifetime::Timeframe); + outputs.emplace_back("ITS3", "VERTICESROF", 0, Lifetime::Timeframe); + outputs.emplace_back("ITS3", "IRFRAMES", 0, Lifetime::Timeframe); + + if (useMC) { + inputs.emplace_back("labels", "ITS3", "CLUSTERSMCTR", 0, Lifetime::Timeframe); + inputs.emplace_back("MC2ROframes", "ITS3", "CLUSTERSMC2ROF", 0, Lifetime::Timeframe); + outputs.emplace_back("ITS3", "VERTICESMCTR", 0, Lifetime::Timeframe); + outputs.emplace_back("ITS3", "TRACKSMCTR", 0, Lifetime::Timeframe); + outputs.emplace_back("ITS3", "ITSTrackMC2ROF", 0, Lifetime::Timeframe); + outputs.emplace_back("ITS3", "VERTICES", 0, Lifetime::Timeframe); + } + + return DataProcessorSpec{ + "its3-tracker", + inputs, + outputs, + AlgorithmSpec{adaptFromTask(ggRequest, useMC, trgType, trModeS, dType)}, + Options{}}; +} + +} // namespace its3 +} // namespace o2 diff --git a/Detectors/Upgrades/IT3/workflow/src/digit-reader-workflow.cxx b/Detectors/Upgrades/ITS3/workflow/src/digit-reader-workflow.cxx similarity index 100% rename from Detectors/Upgrades/IT3/workflow/src/digit-reader-workflow.cxx rename to Detectors/Upgrades/ITS3/workflow/src/digit-reader-workflow.cxx diff --git a/Detectors/Upgrades/IT3/workflow/src/digit-writer-workflow.cxx b/Detectors/Upgrades/ITS3/workflow/src/digit-writer-workflow.cxx similarity index 100% rename from Detectors/Upgrades/IT3/workflow/src/digit-writer-workflow.cxx rename to Detectors/Upgrades/ITS3/workflow/src/digit-writer-workflow.cxx diff --git a/Detectors/Upgrades/IT3/workflow/src/its3-reco-workflow.cxx b/Detectors/Upgrades/ITS3/workflow/src/its3-reco-workflow.cxx similarity index 82% rename from Detectors/Upgrades/IT3/workflow/src/its3-reco-workflow.cxx rename to Detectors/Upgrades/ITS3/workflow/src/its3-reco-workflow.cxx index b202cf67b81f1..5e8c44808e8a0 100644 --- a/Detectors/Upgrades/IT3/workflow/src/its3-reco-workflow.cxx +++ b/Detectors/Upgrades/ITS3/workflow/src/its3-reco-workflow.cxx @@ -44,7 +44,9 @@ void customize(std::vector& workflowOptions) {"clusters-from-upstream", o2::framework::VariantType::Bool, false, {"clusters will be provided from upstream, skip clusterizer"}}, {"disable-root-output", o2::framework::VariantType::Bool, false, {"do not write output root files"}}, {"disable-mc", o2::framework::VariantType::Bool, false, {"disable MC propagation even if available"}}, + {"select-with-triggers", o2::framework::VariantType::String, "none", {"use triggers to prescale processed ROFs: phys, trd, none"}}, {"tracking-mode", o2::framework::VariantType::String, "sync", {"sync,async,cosmics"}}, + {"disable-tracking", o2::framework::VariantType::Bool, false, {"disable tracking step"}}, {"entropy-encoding", o2::framework::VariantType::Bool, false, {"produce entropy encoded data"}}, {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings"}}, {"gpuDevice", o2::framework::VariantType::Int, 1, {"use gpu device: CPU=1,CUDA=2,HIP=3 (default: CPU)"}}}; @@ -62,15 +64,28 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) // Update the (declared) parameters if changed from the command line auto useMC = !configcontext.options().get("disable-mc"); auto trmode = configcontext.options().get("tracking-mode"); + auto selTrig = configcontext.options().get("select-with-triggers"); auto gpuDevice = static_cast(configcontext.options().get("gpuDevice")); auto extDigits = configcontext.options().get("digits-from-upstream"); auto extClusters = configcontext.options().get("clusters-from-upstream"); auto disableRootOutput = configcontext.options().get("disable-root-output"); - auto eencode = configcontext.options().get("entropy-encoding"); + if (configcontext.options().get("disable-tracking")) { + trmode = ""; + } + std::transform(trmode.begin(), trmode.end(), trmode.begin(), [](unsigned char c) { return std::tolower(c); }); o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); - - auto wf = o2::its3::reco_workflow::getWorkflow(useMC, trmode, gpuDevice, extDigits, extClusters, disableRootOutput, eencode); + int trType = 0; + if (!selTrig.empty() && selTrig != "none") { + if (selTrig == "phys") { + trType = 1; + } else if (selTrig == "trd") { + trType = 2; + } else { + LOG(fatal) << "Unknown trigger type requested for events prescaling: " << selTrig; + } + } + auto wf = o2::its3::reco_workflow::getWorkflow(useMC, trmode, gpuDevice, extDigits, extClusters, disableRootOutput, trType); // configure dpl timer to inject correct firstTForbit: start from the 1st orbit of TF containing 1st sampled orbit o2::raw::HBFUtilsInitializer hbfIni(configcontext, wf); diff --git a/Steer/DigitizerWorkflow/CMakeLists.txt b/Steer/DigitizerWorkflow/CMakeLists.txt index 6f9a8d314b7ad..7ac1701c9764e 100644 --- a/Steer/DigitizerWorkflow/CMakeLists.txt +++ b/Steer/DigitizerWorkflow/CMakeLists.txt @@ -19,7 +19,7 @@ o2_add_executable(digitizer-workflow src/GRPUpdaterSpec.cxx src/HMPIDDigitizerSpec.cxx src/ITSMFTDigitizerSpec.cxx - # src/ITS3DigitizerSpec.cxx + src/ITS3DigitizerSpec.cxx src/MCHDigitizerSpec.cxx src/MIDDigitizerSpec.cxx src/PHOSDigitizerSpec.cxx @@ -64,8 +64,8 @@ o2_add_executable(digitizer-workflow O2::ZDCSimulation O2::ZDCWorkflow O2::DetectorsRaw - # O2::ITS3Simulation - # O2::ITS3Workflow + O2::ITS3Simulation + O2::ITS3Workflow ) else() o2_add_executable(digitizer-workflow diff --git a/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.cxx b/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.cxx index 5a8af269ff3a9..9f03bd6733729 100644 --- a/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.cxx +++ b/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.cxx @@ -27,7 +27,7 @@ #include "ITS3Simulation/Digitizer.h" #include "ITSMFTSimulation/DPLDigitizerParam.h" #include "ITSMFTBase/DPLAlpideParam.h" -#include "ITS3Base/GeometryTGeo.h" +#include "ITSBase/GeometryTGeo.h" #include #include #include @@ -81,7 +81,7 @@ class ITS3DPLDigitizerTask : BaseDPLDigitizer << " RO mode"; // configure digitizer - o2::its3::GeometryTGeo* geom = o2::its3::GeometryTGeo::Instance(); + o2::its::GeometryTGeo* geom = o2::its::GeometryTGeo::Instance(); geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); // make sure L2G matrices are loaded mDigitizer.setGeometry(geom); diff --git a/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx b/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx index 24af2c43f8869..65029fbde582a 100644 --- a/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx +++ b/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx @@ -40,11 +40,11 @@ #include "ITSMFTDigitizerSpec.h" #include "ITSMFTWorkflow/DigitWriterSpec.h" -// #ifdef ENABLE_UPGRADES -// // for ITS3 -// #include "ITS3DigitizerSpec.h" -// #include "ITS3Workflow/DigitWriterSpec.h" -// #endif +#ifdef ENABLE_UPGRADES +// for ITS3 +#include "ITS3DigitizerSpec.h" +#include "ITS3Workflow/DigitWriterSpec.h" +#endif // for TOF #include "TOFDigitizerSpec.h" @@ -590,16 +590,16 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) writerSpecs.emplace_back(o2::itsmft::getITSDigitWriterSpec(mctruth)); } - // #ifdef ENABLE_UPGRADES - // // the ITS3 part - // if (isEnabled(o2::detectors::DetID::IT3)) { - // detList.emplace_back(o2::detectors::DetID::IT3); - // // connect the ITS digitization - // specs.emplace_back(o2::its3::getITS3DigitizerSpec(fanoutsize++, mctruth)); - // // // connect ITS digit writer - // specs.emplace_back(o2::its3::getITS3DigitWriterSpec(mctruth)); - // } - // #endif + #ifdef ENABLE_UPGRADES + // the ITS3 part + if (isEnabled(o2::detectors::DetID::IT3)) { + detList.emplace_back(o2::detectors::DetID::IT3); + // connect the ITS digitization + specs.emplace_back(o2::its3::getITS3DigitizerSpec(fanoutsize++, mctruth)); + // // connect ITS digit writer + specs.emplace_back(o2::its3::getITS3DigitWriterSpec(mctruth)); + } + #endif // the MFT part if (isEnabled(o2::detectors::DetID::MFT)) {