diff --git a/Detectors/ITSMFT/common/simulation/include/ITSMFTSimulation/AlpideSimResponse.h b/Detectors/ITSMFT/common/simulation/include/ITSMFTSimulation/AlpideSimResponse.h index 92656a16257a1..5714b51d5aa45 100644 --- a/Detectors/ITSMFT/common/simulation/include/ITSMFTSimulation/AlpideSimResponse.h +++ b/Detectors/ITSMFT/common/simulation/include/ITSMFTSimulation/AlpideSimResponse.h @@ -38,7 +38,7 @@ class AlpideRespSimMat static int constexpr getNPix() { return NPix; } AlpideRespSimMat() = default; - ~AlpideRespSimMat() = default; + virtual ~AlpideRespSimMat() = default; void adopt(const AlpideRespSimMat& src, bool flipRow = false, bool flipCol = false) { @@ -69,7 +69,7 @@ class AlpideRespSimMat private: std::array data; - ClassDefNV(AlpideRespSimMat, 1); + ClassDef(AlpideRespSimMat, 1); }; /* @@ -91,6 +91,7 @@ class AlpideSimResponse int getDepthBin(float pos) const; std::string composeDataName(int colBin, int rowBin); + protected: int mNBinCol = 0; /// number of bins in X(col direction) int mNBinRow = 0; /// number of bins in Y(row direction) int mNBinDpt = 0; /// number of bins in Z(sensor dept) diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h index f8d4a784120a0..fbf9a59e6da4b 100644 --- a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h +++ b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h @@ -12,12 +12,11 @@ /// \file SegmentationMosaix.h /// \brief Definition of the SegmentationMosaix class /// \author felix.schlepper@cern.ch +/// \author chunzheng.wang@cern.ch #ifndef ALICEO2_ITS3_SEGMENTATIONMOSAIX_H_ #define ALICEO2_ITS3_SEGMENTATIONMOSAIX_H_ -#include - #include "MathUtils/Cartesian.h" #include "ITS3Base/SpecsV2.h" @@ -43,24 +42,22 @@ class SegmentationMosaix // 3. The detector coordinate system. Defined by the row and column segmentation // defined at the upper edge in the flat coord. - // row,col=0 - // | - // v - // x----------------------x - // | | | - // | | | - // | | | ^ x - // | | | | - // | | | | - // | | | | - // |-----------X----------| X marks (x,z)=(0,0) X----> z - // | | | + // O----------------------| // | | | + // | | | ^ x + // | | | | + // | | | | + // | | | | + // | | | X----> z X marks (x,z)=(0,0) + // |-----------X----------| + // | | | O----> col O marks (row,col)=(0,0) + // | | | | + // | | | | + // | | | v + // | | | row // | | | - // | | | - // | | | - // | | | - // x----------------------x + // |----------------------| + public: constexpr SegmentationMosaix(int layer) : mRadius(static_cast(constants::radiiMiddle[layer])) {} constexpr ~SegmentationMosaix() = default; @@ -79,7 +76,6 @@ class SegmentationMosaix static constexpr float PitchCol{constants::pixelarray::pixels::mosaix::pitchZ}; static constexpr float PitchRow{constants::pixelarray::pixels::mosaix::pitchX}; static constexpr float SensorLayerThickness{constants::totalThickness}; - static constexpr float NominalYShift{constants::nominalYShift}; /// Transformation from the curved surface to a flat surface. /// Additionally a shift in the flat coordinates must be applied because @@ -102,10 +98,10 @@ class SegmentationMosaix // stack float dist = std::hypot(xCurved, yCurved); float phi = std::atan2(yCurved, xCurved); - xFlat = (mRadius * phi) - WidthH; // the y position is in the silicon volume however we need the chip volume (silicon+metalstack) // this is accounted by a y shift - yFlat = dist - mRadius + NominalYShift; + xFlat = WidthH - mRadius * phi; + yFlat = dist - mRadius; } /// Transformation from the flat surface to a curved surface @@ -122,11 +118,12 @@ class SegmentationMosaix { // MUST align the flat surface with the curved surface with the original pixel array is on and account for metal // stack + float dist = yFlat + mRadius; + float phi = (WidthH - xFlat) / mRadius; // the y position is in the chip volume however we need the silicon volume // this is accounted by a -y shift - float dist = yFlat - NominalYShift + mRadius; - xCurved = dist * std::cos((xFlat + WidthH) / mRadius); - yCurved = dist * std::sin((xFlat + WidthH) / mRadius); + xCurved = dist * std::cos(phi); + yCurved = dist * std::sin(phi); } /// Transformation from Geant detector centered local coordinates (cm) to @@ -142,8 +139,11 @@ class SegmentationMosaix /// \param int iCol Detector z cell coordinate. constexpr bool localToDetector(float const xRow, float const zCol, int& iRow, int& iCol) const noexcept { + if (!isValidLoc(xRow, zCol)) { + return false; + } localToDetectorUnchecked(xRow, zCol, iRow, iCol); - if (!isValid(iRow, iCol)) { + if (!isValidDet(iRow, iCol)) { iRow = iCol = -1; return false; } @@ -167,49 +167,54 @@ class SegmentationMosaix /// center of the sensitive volume. /// If iRow and or iCol is outside of the segmentation range a value of -0.5*Dx() /// or -0.5*Dz() is returned. - constexpr bool detectorToLocal(int const iRow, int const iCol, float& xRow, float& zCol) const noexcept + bool detectorToLocal(float const row, float const col, float& xRow, float& zCol) const noexcept { - if (!isValid(iRow, iCol)) { + if (!isValidDet(row, col)) { return false; } - detectorToLocalUnchecked(iRow, iCol, xRow, zCol); - return isValid(xRow, zCol); + detectorToLocalUnchecked(row, col, xRow, zCol); + return isValidLoc(xRow, zCol); } // Same as detectorToLocal w.o. checks. // We position ourself in the middle of the pixel. - constexpr void detectorToLocalUnchecked(int const iRow, int const iCol, float& xRow, float& zCol) const noexcept + void detectorToLocalUnchecked(float const row, float const col, float& xRow, float& zCol) const noexcept { - xRow = -(static_cast(iRow) + 0.5f) * PitchRow + WidthH; - zCol = (static_cast(iCol) + 0.5f) * PitchCol - LengthH; + xRow = -(row + 0.5f) * PitchRow + WidthH; + zCol = (col + 0.5f) * PitchCol - LengthH; } - bool detectorToLocal(int const row, int const col, math_utils::Point3D& loc) const noexcept + bool detectorToLocal(float const row, float const col, math_utils::Point3D& loc) const noexcept { float xRow{0.}, zCol{0.}; if (!detectorToLocal(row, col, xRow, zCol)) { return false; } - loc.SetCoordinates(xRow, NominalYShift, zCol); + loc.SetCoordinates(xRow, 0.0f, zCol); return true; } - void detectorToLocalUnchecked(int const row, int const col, math_utils::Point3D& loc) const noexcept + void detectorToLocalUnchecked(float const row, float const col, math_utils::Point3D& loc) const noexcept { float xRow{0.}, zCol{0.}; detectorToLocalUnchecked(row, col, xRow, zCol); - loc.SetCoordinates(xRow, NominalYShift, zCol); + loc.SetCoordinates(xRow, 0.0f, zCol); } private: + // Check local coordinates (cm) validity. template - [[nodiscard]] constexpr bool isValid(T const row, T const col) const noexcept + constexpr bool isValidLoc(T const x, T const z) const noexcept { - if constexpr (std::is_floating_point_v) { // compares in local coord. - return (-WidthH < row && row < WidthH && -LengthH < col && col < LengthH); - } else { // compares in rows/cols - return !static_cast(row < 0 || row >= static_cast(NRows) || col < 0 || col >= static_cast(NCols)); - } + return (-WidthH < x && x < WidthH && -LengthH < z && z < LengthH); + } + + // Check detector coordinates validity. + template + constexpr bool isValidDet(T const row, T const col) const noexcept + { + return (row >= 0 && row < static_cast(NRows) && + col >= 0 && col < static_cast(NCols)); } float mRadius; diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h index fedaad9182cce..83db7632e72f4 100644 --- a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h +++ b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h @@ -134,7 +134,6 @@ constexpr std::array radii{19.0006 * mm, 25.228 * mm, 31.4554 * constexpr std::array radiiInner{radii[0] - silicon::thicknessIn, radii[1] - silicon::thicknessIn, radii[2] - silicon::thicknessIn}; // inner silicon radius constexpr std::array radiiOuter{radii[0] + silicon::thicknessOut, radii[1] + silicon::thicknessOut, radii[2] + silicon::thicknessOut}; // outer silicon radius constexpr std::array radiiMiddle{(radiiInner[0] + radiiOuter[0]) / 2., (radiiInner[1] + radiiOuter[1]) / 2., (radiiInner[2] + radiiOuter[2]) / 2.}; // middle silicon radius -constexpr double nominalYShift{-metalstack::thickness / 2.}; // shift to position in silicion volume to the chip volume (silicon+metalstack) // extra information of pixels and their response functions namespace pixelarray::pixels diff --git a/Detectors/Upgrades/ITS3/macros/test/CMakeLists.txt b/Detectors/Upgrades/ITS3/macros/test/CMakeLists.txt index 39e435f0ba2e6..cb6812445283c 100644 --- a/Detectors/Upgrades/ITS3/macros/test/CMakeLists.txt +++ b/Detectors/Upgrades/ITS3/macros/test/CMakeLists.txt @@ -22,7 +22,10 @@ its3_add_macro(CompareClusterSize.C) its3_add_macro(CheckMosaixSegment.C) its3_add_macro(CheckMosaixSegmentTrans.C) its3_add_macro(CompareClustersAndDigits.C) +its3_add_macro(CompareClustersAndDigitsOnChip.C) its3_add_macro(CheckROFs.C) its3_add_macro(CheckTileNumbering.C) its3_add_macro(CreateITS3StaticDeadMap.C) its3_add_macro(TestSensorGeometry.C) +its3_add_macro(CorrTracksClusters.C) +its3_add_macro(CheckChipResponseFile.C) diff --git a/Detectors/Upgrades/ITS3/macros/test/CheckChipResponseFile.C b/Detectors/Upgrades/ITS3/macros/test/CheckChipResponseFile.C new file mode 100644 index 0000000000000..996a99d87ecbc --- /dev/null +++ b/Detectors/Upgrades/ITS3/macros/test/CheckChipResponseFile.C @@ -0,0 +1,192 @@ +// 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 CheckChipResponseFile.C +/// \brief Simple macro to check the chip response files + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#include +#include +#include +#include +#include + +#define ENABLE_UPGRADES +#include "ITSMFTSimulation/AlpideSimResponse.h" + +#include "ITS3Base/SegmentationMosaix.h" +#include "fairlogger/Logger.h" +#endif + +using SegmentationMosaix = o2::its3::SegmentationMosaix; + +double um2cm(double um) { return um * 1e-4; } +double cm2um(double cm) { return cm * 1e+4; } + +o2::itsmft::AlpideSimResponse *mAlpSimResp0 = nullptr, + *mAlpSimResp1 = nullptr, + *mAptSimResp1 = nullptr; + +o2::itsmft::AlpideSimResponse* loadResponse(const std::string& fileName, const std::string& respName) +{ + TFile* f = TFile::Open(fileName.data()); + if (!f) { + std::cerr << fileName << " not found" << std::endl; + return nullptr; + } + auto resp = (o2::itsmft::AlpideSimResponse*)f->Get(respName.data()); + if (!resp) + std::cerr << respName << " not found in " << fileName << std::endl; + return resp; +} + +void LoadRespFunc() +{ + std::string AptsFile = "$(O2_ROOT)/share/Detectors/Upgrades/ITS3/data/ITS3ChipResponseData/APTSResponseData.root"; + std::string AlpideFile = "$(O2_ROOT)/share/Detectors/ITSMFT/data/AlpideResponseData/AlpideResponseData.root"; + + mAlpSimResp0 = loadResponse(AlpideFile, "response0"); // Vbb=0V + LOG(info) << "ALPIDE Vbb=0V response" << std::endl; + mAlpSimResp0->print(); + mAlpSimResp1 = loadResponse(AlpideFile, "response1"); // Vbb=-3V + LOG(info) << "ALPIDE Vbb=-3V response" << std::endl; + mAlpSimResp1->print(); + mAptSimResp1 = loadResponse(AptsFile, "response1"); // APTS + LOG(info) << "APTS response" << std::endl; + mAptSimResp1->print(); +} + +std::vector getCollectionSeediciencies(o2::itsmft::AlpideSimResponse* resp, + const std::vector& depths) +{ + std::vector seed; + bool flipRow = false, flipCol = false; + for (auto depth : depths) { + auto rspmat = resp->getResponse(0.0, 0.0, + um2cm(depth) + resp->getDepthMin() + 1.e-9, + flipRow, flipCol); + seed.push_back(rspmat ? rspmat->getValue(2, 2) : 0.f); + } + return seed; +} + +std::vector getShareValues(o2::itsmft::AlpideSimResponse* resp, + const std::vector& depths) +{ + std::vector share; + bool flipRow = false, flipCol = false; + for (auto depth : depths) { + auto rspmat = resp->getResponse(0.0, 0.0, + um2cm(depth) + resp->getDepthMin() + 1.e-9, + flipRow, flipCol); + float s = 0; + int npix = resp->getNPix(); + if (rspmat) { + for (int i = 0; i < npix; ++i) + for (int j = 0; j < npix; ++j) + if (!(i == npix / 2 && j == npix / 2)) + s += rspmat->getValue(i, j); + } + share.push_back(s); + } + return share; +} + +std::vector getEffValues(o2::itsmft::AlpideSimResponse* resp, + const std::vector& depths) +{ + std::vector all; + bool flipRow = false, flipCol = false; + for (auto depth : depths) { + auto rspmat = resp->getResponse(0.0, 0.0, + um2cm(depth) + resp->getDepthMin() + 1.e-9, + flipRow, flipCol); + float s = 0; + int npix = resp->getNPix(); + if (rspmat) { + for (int i = 0; i < npix; ++i) + for (int j = 0; j < npix; ++j) + s += rspmat->getValue(i, j); + } + all.push_back(s); + } + return all; +} + +void CheckChipResponseFile() +{ + LoadRespFunc(); + LOG(info) << "Response function loaded" << std::endl; + + std::vector vecDepth(50); + for (int i = 0; i < 50; ++i) + vecDepth[i] = i; + + int colors[] = {kOrange + 7, kRed + 1, kAzure + 4}; + struct RespInfo { + o2::itsmft::AlpideSimResponse* resp; + std::string title; + int color; + }; + std::vector responses = { + {mAptSimResp1, "APTS", colors[0]}, + {mAlpSimResp0, "ALPIDE Vbb=0V", colors[1]}, + {mAlpSimResp1, "ALPIDE Vbb=-3V", colors[2]}}; + + TCanvas* c1 = new TCanvas("c1", "c1", 800, 600); + TH1* frame = c1->DrawFrame(-1, -0.049, 50, 1.049); + frame->SetTitle(";Depth(um);Charge Collection Seed / Share / Eff"); + TLegend* leg = new TLegend(0.15, 0.5, 0.4, 0.85); + leg->SetFillStyle(0); + leg->SetBorderSize(0); + + for (auto& r : responses) { + if (!r.resp) + continue; + auto seed = getCollectionSeediciencies(r.resp, vecDepth); + auto shr = getShareValues(r.resp, vecDepth); + auto all = getEffValues(r.resp, vecDepth); + + TGraph* grSeed = new TGraph(vecDepth.size(), vecDepth.data(), seed.data()); + grSeed->SetTitle(Form("%s seed", r.title.c_str())); + grSeed->SetLineColor(r.color); + grSeed->SetLineWidth(2); + grSeed->SetMarkerColor(r.color); + grSeed->SetMarkerStyle(kFullCircle); + grSeed->SetMarkerSize(0.8); + grSeed->Draw("SAME LP"); + leg->AddEntry(grSeed, Form("%s seed", r.title.c_str()), "lp"); + + TGraph* grShare = new TGraph(vecDepth.size(), vecDepth.data(), shr.data()); + grShare->SetLineColor(r.color); + grShare->SetLineWidth(2); + grShare->SetMarkerColor(r.color); + grShare->SetMarkerStyle(kOpenSquare); + grShare->SetMarkerSize(1); + grShare->Draw("SAME LP"); + leg->AddEntry(grShare, Form("%s share", r.title.c_str()), "p"); + + TGraph* grEff = new TGraph(vecDepth.size(), vecDepth.data(), all.data()); + grEff->SetLineColor(r.color); + grEff->SetLineWidth(2); + grEff->SetMarkerColor(r.color); + grEff->SetMarkerStyle(kFullDiamond); + grEff->SetMarkerSize(1); + grEff->Draw("SAME LP"); + leg->AddEntry(grEff, Form("%s eff", r.title.c_str()), "p"); + } + leg->Draw(); + + c1->SaveAs("ChipResponse.pdf"); +} diff --git a/Detectors/Upgrades/ITS3/macros/test/CheckDigitsITS3.C b/Detectors/Upgrades/ITS3/macros/test/CheckDigitsITS3.C index 1dc4a4e2d6b47..240b1bd344af5 100644 --- a/Detectors/Upgrades/ITS3/macros/test/CheckDigitsITS3.C +++ b/Detectors/Upgrades/ITS3/macros/test/CheckDigitsITS3.C @@ -80,8 +80,6 @@ void CheckDigitsITS3(std::string digifile = "it3digits.root", std::string hitfil int nevD = digTree->GetEntries(); // digits in cont. readout may be grouped as few events per entry - int lastReadHitEv = -1; - int nDigitReadIB{0}, nDigitReadOB{0}; int nDigitFilledIB{0}, nDigitFilledOB{0}; diff --git a/Detectors/Upgrades/ITS3/macros/test/CompareClustersAndDigitsOnChip.C b/Detectors/Upgrades/ITS3/macros/test/CompareClustersAndDigitsOnChip.C new file mode 100644 index 0000000000000..310be8c5858ef --- /dev/null +++ b/Detectors/Upgrades/ITS3/macros/test/CompareClustersAndDigitsOnChip.C @@ -0,0 +1,579 @@ +// 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 CompareClustersAndDigitsOnChip.C +/// \brief Macro to compare ITS3 clusters and digits on a pixel array, + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +#define ENABLE_UPGRADES +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsITSMFT/Digit.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "DetectorsCommonDataFormats/DetectorNameConf.h" +#include "ITS3Base/SegmentationMosaix.h" +#include "ITS3Base/SpecsV2.h" +#include "ITS3Reconstruction/TopologyDictionary.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsITSMFT/ClusterTopology.h" +#include "ITSBase/GeometryTGeo.h" +#include "ITSMFTBase/SegmentationAlpide.h" +#include "ITSMFTSimulation/Hit.h" +#include "MathUtils/Cartesian.h" +#include "MathUtils/Utils.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "SimulationDataFormat/ConstMCTruthContainer.h" +#include "SimulationDataFormat/IOMCTruthContainerView.h" + +struct Data { + TH2F* pixelArray; + TGraph* hitS; + TGraph* hitM; + TGraph* hitE; + TGraph* clusS; + TGraph* cog; + TLegend* leg; + std::vector* vClusBox; + void clear() + { + delete pixelArray; + delete hitS; + delete hitM; + delete hitE; + delete clusS; + delete cog; + delete leg; + for (auto& b : *vClusBox) { + delete b; + } + delete vClusBox; + } +}; + +void CompareClustersAndDigitsOnChip(std::string clusfile = "o2clus_its.root", + std::string digifile = "it3digits.root", + std::string dictfile = "", + std::string hitfile = "o2sim_HitsIT3.root", + std::string inputGeom = "o2sim_geometry.root", + bool batch = true) +{ + TH1::AddDirectory(kFALSE); + gROOT->SetBatch(batch); + gStyle->SetPalette(kRainBow); + gStyle->SetOptStat(0); + + using namespace o2::base; + using namespace o2::its; + using o2::itsmft::Hit; + using Segmentation = o2::itsmft::SegmentationAlpide; + using o2::itsmft::ClusterTopology; + using o2::itsmft::CompClusterExt; + 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; + + std::array mMosaixSegmentations{0, 1, 2}; + + // 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 + const int nChips = gman->getNumberOfChips(); + + LOGP(info, "Total number of chips is {} in ITS3 (IB and OB)", nChips); + + // Create all plots + LOGP(info, "Selecting chips to be visualised"); + std::set selectedChips; + std::map> chipGroups; + + for (int chipID{0}; chipID < nChips; ++chipID) { + TString tpath = gman->getMatrixPath(chipID); + std::string path = tpath.Data(); + + std::vector tokens; + std::istringstream iss(path); + std::string token; + while (std::getline(iss, token, '/')) { + if (!token.empty()) { + tokens.push_back(token); + } + } + + std::string segmentName, staveName, carbonFormName; + for (const auto& t : tokens) { + if (t.find("ITS3Segment") != std::string::npos) + segmentName = t; + if (t.find("ITSUStave") != std::string::npos) + staveName = t; + if (t.find("ITS3CarbonForm") != std::string::npos) + carbonFormName = t; + } + + std::string groupKey; + if (!segmentName.empty()) { + groupKey = segmentName + "_" + carbonFormName; + } else if (!staveName.empty()) { + groupKey = staveName; + } else { + continue; + } + + chipGroups[groupKey].push_back(chipID); + } + + LOGP(info, "From each IB Segment or OB Stave, 10 chipIDs are uniformly selected"); + LOGP(info, "Selected chipID: "); + for (auto& [groupName, ids] : chipGroups) { + std::vector sampled; + if (ids.size() <= 10) { + for (auto id : ids) { + selectedChips.insert(id); + sampled.push_back(id); + } + } else { + for (int i{0}; i < 10; ++i) { + int idx = i * (ids.size() - 1) / 9; // 9 intervals for 10 points + int id = ids[idx]; + if (selectedChips.insert(id).second) { + sampled.push_back(id); + } + } + } + + std::ostringstream oss; + std::string topOrBot = "N/A"; + std::smatch match; + std::regex rgxSegment(R"(Segment(\d+)_(\d+)_ITS3CarbonForm\d+_(\d+))"); + std::regex rgxStave(R"(Stave(\d+)_(\d+))"); + if (std::regex_search(groupName, match, rgxSegment)) { + int layer = std::stoi(match[1]); + int segment = std::stoi(match[2]); + int carbonForm = std::stoi(match[3]); + topOrBot = (carbonForm == 0 ? "TOP" : "BOT"); + oss << topOrBot << " segment " << segment << " at layer " << layer << ": "; + } else if (std::regex_search(groupName, match, rgxStave)) { + int layer = std::stoi(match[1]); + int stave = std::stoi(match[2]); + oss << "Stave " << stave << " at layer " << layer << ": "; + } else { + LOGP(error, "Cannot select the correct chipID in OB or IB"); + return; + } + for (auto id : sampled) { + oss << id << " "; + } + LOG(info) << oss.str(); + } + LOGP(info, "{} selected chips will be visualized and analyzed.", chipGroups.size()); + + // Hits + TFile fileH(hitfile.data()); + auto* hitTree = dynamic_cast(fileH.Get("o2sim")); + std::vector* hitArray = nullptr; + hitTree->SetBranchAddress("IT3Hit", &hitArray); + mc2hitVec.resize(hitTree->GetEntries()); + hitVecPool.resize(hitTree->GetEntries(), nullptr); + + // Digits + TFile* digFile = TFile::Open(digifile.data()); + TTree* digTree = (TTree*)digFile->Get("o2sim"); + std::vector* digArr = nullptr; + digTree->SetBranchAddress("IT3Digit", &digArr); + o2::dataformats::IOMCTruthContainerView* plabels = nullptr; + digTree->SetBranchAddress("IT3DigitMCTruth", &plabels); + + // Clusters + TFile fileC(clusfile.data()); + auto* clusTree = dynamic_cast(fileC.Get("o2sim")); + std::vector* clusArr = nullptr; + clusTree->SetBranchAddress("ITSClusterComp", &clusArr); + std::vector* patternsPtr = nullptr; + auto pattBranch = clusTree->GetBranch("ITSClusterPatt"); + if (pattBranch != nullptr) { + pattBranch->SetAddress(&patternsPtr); + } + + // Topology dictionary + o2::its3::TopologyDictionary dict; + bool hasAvailableDict = false; + if (!dictfile.empty()) { + std::ifstream file(dictfile.c_str()); + if (file.good()) { + LOGP(info, "Running with external topology dictionary: {}", dictfile); + dict.readFromFile(dictfile); + LOGP(info, "The IB dictionary size is {}, and the OB dictionary size is {}", dict.getSize(true), dict.getSize(false)); + hasAvailableDict = dict.getSize(true) != 0 && dict.getSize(false) != 0; + if (hasAvailableDict) { + LOGP(info, "Dictionaries is vaild."); + } else { + LOGP(info, "Dictionaries is NOT vaild!"); + } + } else { + LOGP(info, "Cannot open dictionary file: {}. Running without external dictionary!", dictfile); + dictfile = ""; + } + } else { + LOGP(info, "Running without external topology dictionary!"); + } + + // ROFrecords + std::vector rofRecVec, *rofRecVecP = &rofRecVec; + clusTree->SetBranchAddress("ITSClustersROF", &rofRecVecP); + + // Cluster MC labels + o2::dataformats::MCTruthContainer* clusLabArr = nullptr; + std::vector mc2rofVec, *mc2rofVecP = &mc2rofVec; + if ((hitTree != nullptr) && (clusTree->GetBranch("ITSClusterMCTruth") != nullptr)) { + clusTree->SetBranchAddress("ITSClusterMCTruth", &clusLabArr); + clusTree->SetBranchAddress("ITSClustersMC2ROF", &mc2rofVecP); + } + + clusTree->GetEntry(0); + unsigned 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]; + if (mc2rof.rofRecordID < 0) { + continue; // this MC event did not contribute to any ROF + } + for (unsigned int irfd = mc2rof.maxROF - mc2rof.minROF + 1; irfd--;) { + unsigned int irof = mc2rof.rofRecordID + irfd; + if (irof >= nROFRec) { + LOGP(error, "ROF = {} from MC2ROF record is >= N ROFs = {}", irof, nROFRec); + } + if (mcEvMin[irof] > imc) { + mcEvMin[irof] = imc; + } + if (mcEvMax[irof] < imc) { + mcEvMax[irof] = imc; + } + } + } + + // Create all plots + LOGP(info, "Creating plots"); + std::unordered_map data; + auto initData = [&](int chipID, Data& dat) { + if (dat.pixelArray) + return; + + int nCol{0}, nRow{0}; + float lengthPixArr{0}, widthPixArr{0}; + bool isIB = o2::its3::constants::detID::isDetITS3(chipID); + int layer = gman->getLayer(chipID); + if (isIB) { + nCol = o2::its3::SegmentationMosaix::NCols; + nRow = o2::its3::SegmentationMosaix::NRows; + lengthPixArr = o2::its3::constants::pixelarray::pixels::mosaix::pitchZ * nCol; + widthPixArr = o2::its3::constants::pixelarray::pixels::mosaix::pitchX * nRow; + } else { + nCol = o2::itsmft::SegmentationAlpide::NCols; + nRow = o2::itsmft::SegmentationAlpide::NRows; + lengthPixArr = o2::itsmft::SegmentationAlpide::PitchCol * nCol; + widthPixArr = o2::itsmft::SegmentationAlpide::PitchRow * nRow; + } + + dat.pixelArray = new TH2F(Form("histSensor_%d", chipID), Form("SensorID=%d;z(cm);x(cm)", chipID), + nCol, -0.5 * lengthPixArr, 0.5 * lengthPixArr, + nRow, -0.5 * widthPixArr, 0.5 * widthPixArr); + dat.hitS = new TGraph(); + dat.hitS->SetMarkerStyle(kFullTriangleDown); + dat.hitS->SetMarkerColor(kGreen); + dat.hitM = new TGraph(); + dat.hitM->SetMarkerStyle(kFullCircle); + dat.hitM->SetMarkerColor(kGreen + 3); + dat.hitE = new TGraph(); + dat.hitE->SetMarkerStyle(kFullTriangleUp); + dat.hitE->SetMarkerColor(kGreen + 5); + dat.clusS = new TGraph(); + dat.clusS->SetMarkerStyle(kFullSquare); + dat.clusS->SetMarkerColor(kBlue); + dat.cog = new TGraph(); + dat.cog->SetMarkerStyle(kFullDiamond); + dat.cog->SetMarkerColor(kRed); + dat.leg = new TLegend(0.7, 0.7, 0.92, 0.92); + dat.leg->AddEntry(dat.hitS, "Hit Start"); + dat.leg->AddEntry(dat.hitM, "Hit Middle"); + dat.leg->AddEntry(dat.hitE, "Hit End"); + dat.leg->AddEntry(dat.clusS, "Cluster Start"); + dat.leg->AddEntry(dat.cog, "Cluster COG"); + dat.vClusBox = new std::vector; + }; + + LOGP(info, "Filling digits"); + for (int iDigit{0}; digTree->LoadTree(iDigit) >= 0; ++iDigit) { + digTree->GetEntry(iDigit); + for (const auto& digit : *digArr) { + const auto chipID = digit.getChipIndex(); + if (!selectedChips.count(chipID)) + continue; + const auto layer = gman->getLayer(chipID); + bool isIB = layer < 3; + float locDigiX{0}, locDigiZ{0}; + if (isIB) { + mMosaixSegmentations[layer].detectorToLocal(digit.getRow(), digit.getColumn(), locDigiX, locDigiZ); + } else { + o2::itsmft::SegmentationAlpide::detectorToLocal(digit.getRow(), digit.getColumn(), locDigiX, locDigiZ); + } + auto& dat = data[chipID]; + initData(chipID, dat); + data[chipID].pixelArray->Fill(locDigiZ, locDigiX); + } + } + + LOGP(info, "Building min and max MC events used by each ROF, total ROFs {}", nROFRec); + auto pattIt = patternsPtr->cbegin(); + bool isAllPattIDInvaild{true}; + for (unsigned int irof{0}; irof < nROFRec; irof++) { + const auto& rofRec = rofRecVec[irof]; + // >> read and map MC events contributing to this ROF + for (int im = mcEvMin[irof]; im <= mcEvMax[irof]; im++) { + if (hitVecPool[im] == nullptr) { + 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); + } + } + } + + // Clusters in 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]; + const auto chipID = cluster.getSensorID(); + if (!selectedChips.count(chipID)) { + // Even if not selected, advance pattIt if patternID is InvalidPatternID + if (cluster.getPatternID() == o2::itsmft::CompCluster::InvalidPatternID) { + o2::itsmft::ClusterPattern::skipPattern(pattIt); + } + continue; + } + const auto pattID = cluster.getPatternID(); + const bool isIB = o2::its3::constants::detID::isDetITS3(chipID); + const auto layer = gman->getLayer(chipID); + auto& dat = data[chipID]; + initData(chipID, dat); + o2::itsmft::ClusterPattern pattern; + // Pattern extraction + if (cluster.getPatternID() != o2::itsmft::CompCluster::InvalidPatternID) { + isAllPattIDInvaild = false; + if (!hasAvailableDict) { + LOGP(error, "Encountered pattern ID {}, which is not equal to the invalid pattern ID {}", cluster.getPatternID(), o2::itsmft::CompCluster::InvalidPatternID); + LOGP(error, "Clusters have already been generated with a dictionary which was not provided properly!"); + return; + } + if (dict.isGroup(cluster.getPatternID(), isIB)) { + pattern.acquirePattern(pattIt); + } else { + pattern = dict.getPattern(cluster.getPatternID(), isIB); + } + } else { + pattern.acquirePattern(pattIt); + } + + // Hits + const auto& lab = (clusLabArr->getLabels(clEntry))[0]; + if (!lab.isValid()) + continue; + const 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()) + continue; + o2::math_utils::Point3D locHMiddle; + const auto& hit = (*hitArray)[hitEntry->second]; + auto locHEnd = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); + auto locHStart = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart()); + if (isIB) { + float xFlat{0.}, yFlat{0.}; + mMosaixSegmentations[layer].curvedToFlat(locHEnd.X(), locHEnd.Y(), xFlat, yFlat); + locHEnd.SetXYZ(xFlat, yFlat, locHEnd.Z()); + mMosaixSegmentations[layer].curvedToFlat(locHStart.X(), locHStart.Y(), xFlat, yFlat); + locHStart.SetXYZ(xFlat, yFlat, locHStart.Z()); + } + locHMiddle.SetXYZ(0.5f * (locHEnd.X() + locHStart.X()), + 0.5f * (locHEnd.Y() + locHStart.Y()), + 0.5f * (locHEnd.Z() + locHStart.Z())); + data[chipID].hitS->AddPoint(locHStart.Z(), locHStart.X()); + data[chipID].hitM->AddPoint(locHMiddle.Z(), locHMiddle.X()); + data[chipID].hitE->AddPoint(locHEnd.Z(), locHEnd.X()); + + // Cluster Start point + float locCluX{0}, locCluZ{0}; + if (isIB) { + mMosaixSegmentations[layer].detectorToLocal(cluster.getRow(), cluster.getCol(), locCluX, locCluZ); + } else { + o2::itsmft::SegmentationAlpide::detectorToLocal(cluster.getRow(), cluster.getCol(), locCluX, locCluZ); + } + data[chipID].clusS->AddPoint(locCluZ, locCluX); + + // COG + o2::math_utils::Point3D locCOG; + // Cluster COG using dictionary (if available) + if (hasAvailableDict && (pattID != o2::itsmft::CompCluster::InvalidPatternID && !dict.isGroup(pattID, isIB))) { + locCOG = dict.getClusterCoordinates(cluster); + } else { + if (isIB) { + locCOG = o2::its3::TopologyDictionary::getClusterCoordinates(cluster, pattern, false); + } else { + locCOG = o2::itsmft::TopologyDictionary::getClusterCoordinates(cluster, pattern, false); + } + } + if (isIB) { + float flatX{0}, flatY{0}; + mMosaixSegmentations[layer].curvedToFlat(locCOG.X(), locCOG.Y(), flatX, flatY); + locCOG.SetCoordinates(flatX, flatY, locCOG.Z()); + } + data[chipID].cog->AddPoint(locCOG.Z(), locCOG.X()); + + // Cluster Box using dictionary if available, otherwise use raw pattern + float lowLeftX{0}, lowLeftZ{0}, topRightX{0}, topRightZ{0}; + // Use dictionary-based cluster box + if (isIB) { + mMosaixSegmentations[layer].detectorToLocal(cluster.getRow(), cluster.getCol(), lowLeftX, lowLeftZ); + mMosaixSegmentations[layer].detectorToLocal(cluster.getRow() + pattern.getRowSpan() - 1, + cluster.getCol() + pattern.getColumnSpan() - 1, + topRightX, topRightZ); + lowLeftX += 0.5 * o2::its3::constants::pixelarray::pixels::mosaix::pitchX; + lowLeftZ -= 0.5 * o2::its3::constants::pixelarray::pixels::mosaix::pitchZ; + topRightX -= 0.5 * o2::its3::constants::pixelarray::pixels::mosaix::pitchX; + topRightZ += 0.5 * o2::its3::constants::pixelarray::pixels::mosaix::pitchZ; + } else { + o2::itsmft::SegmentationAlpide::detectorToLocal(cluster.getRow(), cluster.getCol(), lowLeftX, lowLeftZ); + o2::itsmft::SegmentationAlpide::detectorToLocal(cluster.getRow() + pattern.getRowSpan() - 1, + cluster.getCol() + pattern.getColumnSpan() - 1, + topRightX, topRightZ); + lowLeftX += 0.5 * o2::itsmft::SegmentationAlpide::PitchRow; + lowLeftZ -= 0.5 * o2::itsmft::SegmentationAlpide::PitchCol; + topRightX -= 0.5 * o2::itsmft::SegmentationAlpide::PitchRow; + topRightZ += 0.5 * o2::itsmft::SegmentationAlpide::PitchCol; + } + auto clusBox = new TBox(lowLeftZ, lowLeftX, topRightZ, topRightX); + clusBox->SetFillColorAlpha(0, 0); + clusBox->SetFillStyle(0); + clusBox->SetLineWidth(4); + clusBox->SetLineColor(kBlack); + data[chipID].vClusBox->push_back(clusBox); + } + } + + if (isAllPattIDInvaild) { + LOGP(info, "Verified input cluster file was generated w/o topology dictionary"); + if (!dictfile.empty()) { + LOGP(error, "Non-dictionary cluster file processed by external dictionary! Please adjust input."); + return; + } + } + + LOGP(info, "Writing to root file"); + double x1, y1, x2, y2; + auto oFileOut = TFile::Open("CompareClustersAndDigitsOnChip.root", "RECREATE"); + oFileOut->cd(); + for (int chipID{0}; chipID < nChips; chipID++) { + if (!selectedChips.count(chipID)) + continue; + auto& dat = data[chipID]; + TString tpath = gman->getMatrixPath(chipID); + const std::string cpath{tpath.Data() + 39, tpath.Data() + tpath.Length()}; + const std::filesystem::path p{cpath}; + std::string nestedDir = p.parent_path().string(); + TDirectory* currentDir = oFileOut; + std::istringstream iss(nestedDir); + std::string token; + while (std::getline(iss, token, '/')) { + if (token.empty()) + continue; + TDirectory* nextDir = currentDir->GetDirectory(token.c_str()); + if (!nextDir) { + nextDir = currentDir->mkdir(token.c_str()); + } + if (!nextDir) { + LOGP(error, "Cannot create subdirectory: %s", token.c_str()); + break; + } + currentDir = nextDir; + currentDir->cd(); + } + if (!currentDir) { + LOGP(error, "Failed to create nested directory for chip %d", chipID); + continue; + } + + auto canv = new TCanvas(Form("%s_%d", p.filename().c_str(), chipID)); + canv->SetTitle(Form("%s_%d", p.filename().c_str(), chipID)); + canv->cd(); + gPad->SetGrid(1, 1); + dat.pixelArray->Draw("colz"); + dat.hitS->Draw("p;same"); + dat.hitM->Draw("p;same"); + dat.hitE->Draw("p;same"); + auto arr = new TArrow(); + arr->SetArrowSize(0.01); + for (int i{0}; i < dat.hitS->GetN(); ++i) { + dat.hitS->GetPoint(i, x1, y1); + dat.hitE->GetPoint(i, x2, y2); + arr->DrawArrow(x1, y1, x2, y2); + } + dat.clusS->Draw("p;same"); + if (dat.cog->GetN() != 0) + dat.cog->Draw("p;same"); + for (const auto& clusBox : *dat.vClusBox) { + clusBox->Draw(); + } + dat.leg->Draw(); + canv->SetEditable(false); + + currentDir->WriteTObject(canv, canv->GetName()); + dat.clear(); + delete canv; + delete arr; + printf("\rWriting chip %05d", chipID); + } + printf("\n"); + oFileOut->Write(); + oFileOut->Close(); + LOGP(info, "Finished writing selected chip visualizations."); +} \ No newline at end of file diff --git a/Detectors/Upgrades/ITS3/macros/test/CorrTracksClusters.C b/Detectors/Upgrades/ITS3/macros/test/CorrTracksClusters.C new file mode 100644 index 0000000000000..634d761366920 --- /dev/null +++ b/Detectors/Upgrades/ITS3/macros/test/CorrTracksClusters.C @@ -0,0 +1,638 @@ +// 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 "TEfficiency.h" +#include +#include +#include + +#include "ITSMFTSimulation/Hit.h" +#include "DataFormatsITS/TrackITS.h" +#include "DetectorsBase/Propagator.h" +#include "Field/MagneticField.h" +#include "ITSBase/GeometryTGeo.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCEventHeader.h" +#include "SimulationDataFormat/MCTrack.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "SimulationDataFormat/TrackReference.h" +#include "ITS3Reconstruction/TopologyDictionary.h" +#include "ITSMFTBase/SegmentationAlpide.h" +#include "ITS3Base/SegmentationMosaix.h" + +#include +#include +#include +#include +#endif + +using namespace std; +using namespace o2::itsmft; +using namespace o2::its; +using SegmentationIB = o2::its3::SegmentationMosaix; +using SegmentationOB = o2::itsmft::SegmentationAlpide; +static constexpr int kNLayer = 7; +static constexpr int INVALID_INT = -99; +static constexpr float INVALID_FLOAT = -99.f; + +//______________________________________________________________________________ +// ParticleInfo structure +struct ParticleInfo { + int event{}; + int pdg{}; + float pt{}; + float recpt{}; + float eta{}; + float phi{}; + float pvx{}; + float pvy{}; + float pvz{}; + float dcaxy{}; + float dcaz{}; + int mother{}; + int first{}; + unsigned short clusters = 0u; + unsigned char isReco = 0u; + unsigned char isFake = 0u; + bool isPrimary = false; + unsigned char storedStatus = 2; /// not stored = 2, fake = 1, good = 0 + std::array clusterSize; + std::array clusterPattern; + std::array clusterLocX; + std::array clusterLocZ; + std::array hitLocX; + std::array hitLocY; + std::array hitLocZ; + o2::its::TrackITS track; + ParticleInfo() + { + clusterSize.fill(INVALID_INT); + clusterPattern.fill(INVALID_INT); + clusterLocX.fill(INVALID_FLOAT); + clusterLocZ.fill(INVALID_FLOAT); + hitLocX.fill(INVALID_FLOAT); + hitLocY.fill(INVALID_FLOAT); + hitLocZ.fill(INVALID_FLOAT); + } +}; + +//______________________________________________________________________________ +// Convert curved local coordinates to flat coordinates +void CurvedLocalToFlat(o2::math_utils::Point3D& point, const SegmentationIB& seg) +{ + float xFlat = 0.f, yFlat = 0.f; + seg.curvedToFlat(point.X(), point.Y(), xFlat, yFlat); + point.SetXYZ(xFlat, yFlat, point.Z()); +} + +//______________________________________________________________________________ +// Resolve pattern from patternID and iterator +bool resolvePattern(const o2::itsmft::CompClusterExt& cluster, + decltype(std::declval>().cbegin())& pattIt, + const o2::its3::TopologyDictionary& dict, + bool isIB, + o2::itsmft::ClusterPattern& pattOut) +{ + auto pattID = cluster.getPatternID(); + if (pattID != o2::itsmft::CompCluster::InvalidPatternID) { + if (!dict.getSize(true) && !dict.getSize(false)) { + LOGP(error, "Encountered non-invalid pattern ID {} but dictionary is missing!", pattID); + return false; + } + if (dict.isGroup(pattID, isIB)) { + pattOut.acquirePattern(pattIt); + } else { + pattOut = dict.getPattern(pattID, isIB); + } + } else { + pattOut.acquirePattern(pattIt); + } + return true; +} + +//______________________________________________________________________________ +// Function to analyze reconstructed tracks +void analyzeRecoTracks(TTree* recTree, + const std::vector* recArr, + const std::vector* trkLabArr, + std::vector>& info, + float bz, + ULong_t& unaccounted, + ULong_t& good, + ULong_t& fakes, + ULong_t& total) +{ + unaccounted = good = fakes = total = 0; + for (int frame = 0; frame < recTree->GetEntriesFast(); frame++) { // reco tracks frames + if (recTree->GetEvent(frame) == 0) + continue; + total += trkLabArr->size(); + for (unsigned int iTrack = 0; iTrack < trkLabArr->size(); ++iTrack) { + auto lab = trkLabArr->at(iTrack); + if (!lab.isSet()) { + unaccounted++; + continue; + } + int trackID, evID, srcID; + bool fake; + lab.get(trackID, evID, srcID, fake); + if (evID < 0 || evID >= (int)info.size()) { + unaccounted++; + continue; + } + if (trackID < 0 || trackID >= (int)info[evID].size()) { + unaccounted++; + continue; + } + info[evID][trackID].isReco += !fake; + info[evID][trackID].isFake += fake; + if (recArr->at(iTrack).isBetter(info[evID][trackID].track, 1.e9)) { + info[evID][trackID].storedStatus = fake; + info[evID][trackID].track = recArr->at(iTrack); + float ip[2]{0., 0.}; + info[evID][trackID].track.getImpactParams(info[evID][trackID].pvx, + info[evID][trackID].pvy, + info[evID][trackID].pvz, bz, ip); + info[evID][trackID].dcaxy = ip[0]; + info[evID][trackID].dcaz = ip[1]; + info[evID][trackID].recpt = info[evID][trackID].track.getPt(); + } + fakes += static_cast(fake); + good += static_cast(!fake); + } + } + LOGP(info, "** Some statistics:"); + LOGP(info, "\t- Total number of tracks: {}", total); + LOGP(info, "\t- Total number of tracks not corresponding to particles: {} ({:.2f}%)", unaccounted, unaccounted * 100. / total); + LOGP(info, "\t- Total number of fakes: {} ({:.2f}%)", fakes, fakes * 100. / total); + LOGP(info, "\t- Total number of good: {} ({:.2f}%)", good, good * 100. / total); +} + +//______________________________________________________________________________ +// Read and map hit information from hitTree +void mapHitsForMCEvents(TTree* hitTree, + std::vector*>& hitVecPool, + std::vector>& mc2hitVec, + const std::vector& mcEvMin, + const std::vector& mcEvMax, + size_t nROFRec) +{ + for (unsigned int irof = 0; irof < nROFRec; irof++) { + 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); + } + } + } + } +} + +//______________________________________________________________________________ +// Load geometry and magnetic field information +void loadGeometryAndField(const std::string& magfile, const std::string& inputGeom, float& bz, o2::its::GeometryTGeo*& gman) +{ + o2::base::Propagator::initFieldFromGRP(magfile); + bz = o2::base::Propagator::Instance()->getNominalBz(); + o2::base::GeometryManager::loadGeometry(inputGeom); + 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)); +} + +//______________________________________________________________________________ +// Load topology dictionary +void loadTopologyDictionary(const std::string& dictfile, o2::its3::TopologyDictionary& dict) +{ + std::ifstream iofile(dictfile); + if (iofile.good()) { + LOG(info) << "Running with dictionary: " << dictfile; + dict.readFromFile(dictfile); + } else { + LOG(info) << "Dictionary file not found: " << dictfile; + } +} + +//______________________________________________________________________________ +// Build ROF +void buildMcEvRangePerROF(const std::vector& mc2rofVec, + size_t nROFRec, + std::vector& mcEvMin, + std::vector& mcEvMax) +{ + for (size_t imc = 0; imc < mc2rofVec.size(); ++imc) { + const auto& mc2rof = mc2rofVec[imc]; + if (mc2rof.rofRecordID < 0) + continue; + for (size_t i = mc2rof.minROF; i <= mc2rof.maxROF; ++i) { + if (i >= nROFRec) + continue; + mcEvMin[i] = std::min(mcEvMin[i], static_cast(imc)); + mcEvMax[i] = std::max(mcEvMax[i], static_cast(imc)); + } + } +} + +//______________________________________________________________________________ +// Load Hits data +void prepareHitAccess(const std::string& hitfile, + TTree*& hitTree, + std::vector*>& hitVecPool, + std::vector>& mc2hitVec) +{ + TFile* fHit = TFile::Open(hitfile.data()); + hitTree = (TTree*)fHit->Get("o2sim"); + mc2hitVec.resize(hitTree->GetEntries()); + hitVecPool.resize(hitTree->GetEntries(), nullptr); +} + +void loadCluster(const std::string& clusfile, + TTree*& clusTree, + std::vector*& clusArr, + o2::dataformats::MCTruthContainer*& clusLabArr, + std::vector& mc2rofVec, + std::vector*& patternsPtr, + std::vector& rofRecVec) +{ + // Open file and let it persist + TFile* fileC = TFile::Open(clusfile.data()); + // Get tree + clusTree = dynamic_cast(fileC->Get("o2sim")); + // Cluster array + clusArr = nullptr; + clusTree->SetBranchAddress("ITSClusterComp", &clusArr); + // MC truth + clusLabArr = nullptr; + clusTree->SetBranchAddress("ITSClusterMCTruth", &clusLabArr); + clusTree->SetBranchAddress("ITSClusterPatt", &patternsPtr); + // ROF records + std::vector* rofRecVecP = &rofRecVec; + clusTree->SetBranchAddress("ITSClustersROF", &rofRecVecP); + // MC2ROF mapping + std::vector* mc2rofVecP = &mc2rofVec; + clusTree->SetBranchAddress("ITSClustersMC2ROF", &mc2rofVecP); + clusTree->GetEntry(0); + // After setting all branch addresses, trigger preload of the first entr +} + +//______________________________________________________________________________ +// Load Reconstructed Tracks data +void loadRecoTracks(const std::string& tracfile, + TTree*& recTree, + std::vector*& recArr, + std::vector*& trkLabArr) +{ + TFile* fTrk = TFile::Open(tracfile.data()); + recTree = (TTree*)fTrk->Get("o2sim"); + recTree->SetBranchAddress("ITSTrack", &recArr); + recTree->SetBranchAddress("ITSTrackMCTruth", &trkLabArr); +} + +//______________________________________________________________________________ +// Load MC Track information +void loadMCTrackInfo(const std::string& kinefile, + std::vector>& info, + std::vector*& mcArr, + o2::dataformats::MCEventHeader*& mcEvent, + TTree*& mcTree) +{ + TFile* kineFile = TFile::Open(kinefile.data()); + mcTree = (TTree*)kineFile->Get("o2sim"); + mcTree->SetBranchStatus("*", 0); + mcTree->SetBranchStatus("MCTrack*", 1); + mcTree->SetBranchStatus("MCEventHeader*", 1); + mcTree->SetBranchAddress("MCTrack", &mcArr); + mcTree->SetBranchAddress("MCEventHeader.", &mcEvent); + + int nev = mcTree->GetEntriesFast(); + info.resize(nev); + for (int n = 0; n < nev; n++) { + mcTree->GetEvent(n); + info[n].resize(mcArr->size()); + for (unsigned int mcI = 0; mcI < mcArr->size(); ++mcI) { + auto part = mcArr->at(mcI); + info[n][mcI].pvx = mcEvent->GetX(); + info[n][mcI].pvy = mcEvent->GetY(); + info[n][mcI].pvz = mcEvent->GetZ(); + info[n][mcI].event = n; + info[n][mcI].pdg = part.GetPdgCode(); + info[n][mcI].pt = part.GetPt(); + info[n][mcI].phi = part.GetPhi(); + info[n][mcI].eta = part.GetEta(); + info[n][mcI].isPrimary = part.isPrimary(); + } + } +} + +//______________________________________________________________________________ +// Main function CorrTracksClusters +void CorrTracksClusters(const std::string& tracfile = "o2trac_its.root", + const std::string& clusfile = "o2clus_its.root", + const std::string& kinefile = "o2sim_Kine.root", + const std::string& magfile = "o2sim_grp.root", + const std::string& hitfile = "o2sim_HitsIT3.root", + const std::string& dictfile = "IT3dictionary.root", + const std::string& inputGeom = "", + bool batch = false) +{ + gROOT->SetBatch(batch); + + // Geo and Field + LOGP(info, "Geo and Field loading"); + float bz{0.f}; + o2::its::GeometryTGeo* gman = nullptr; + loadGeometryAndField(magfile, inputGeom, bz, gman); + LOGP(info, "Finished Geo and Field loading"); + + // MC tracks + LOGP(info, "MC Track Info loading"); + std::vector* mcArr = nullptr; + o2::dataformats::MCEventHeader* mcEvent = nullptr; + TTree* mcTree = nullptr; + std::vector> info; + loadMCTrackInfo(kinefile, info, mcArr, mcEvent, mcTree); + LOGP(info, "Finished MC Track Info loading"); + + // Reconstructed tracks + LOGP(info, "Reco Tracks loading"); + TTree* recTree = nullptr; + std::vector* recArr = nullptr; + std::vector* trkLabArr = nullptr; + loadRecoTracks(tracfile, recTree, recArr, trkLabArr); + LOGP(info, "Finished Reco Tracks loading"); + + // Run analyzeRecoTracks + LOGP(info, "Track analysis (analyzeRecoTracks)"); + ULong_t unaccounted{0}, good{0}, fakes{0}, total{0}; + analyzeRecoTracks(recTree, recArr, trkLabArr, info, bz, unaccounted, good, fakes, total); + LOGP(info, "Finished track analysis (analyzeRecoTracks)"); + + // Topology dictionary + LOGP(info, "Topology Dictionary loading"); + o2::its3::TopologyDictionary dict; + loadTopologyDictionary(dictfile, dict); + LOGP(info, "Finished Topology Dictionary loading"); + + // Clusters + LOGP(info, "Cluster Data loading"); + TTree* clusTree = nullptr; + std::vector* clusArr = nullptr; + o2::dataformats::MCTruthContainer* clusLabArr = nullptr; + std::vector* patternsPtr = nullptr; + std::vector mc2rofVec; + std::vector rofRecVec; + loadCluster(clusfile, clusTree, clusArr, clusLabArr, mc2rofVec, patternsPtr, rofRecVec); + LOGP(info, "Finished Cluster Data loading"); + // clusTree->GetEntry(0); + + // Hits + LOGP(info, "Hits loading"); + TTree* hitTree = nullptr; + std::vector*> hitVecPool; + std::vector> mc2hitVec; + prepareHitAccess(hitfile, hitTree, hitVecPool, mc2hitVec); + LOGP(info, "Finished Hits loading"); + + // Build min and max MC events used by each ROF + LOGP(info, "Building MC event ranges"); + std::vector mcEvMin, mcEvMax; + mcEvMin.assign(rofRecVec.size(), hitTree->GetEntries()); + mcEvMax.assign(rofRecVec.size(), -1); + buildMcEvRangePerROF(mc2rofVec, rofRecVec.size(), mcEvMin, mcEvMax); + LOGP(info, "Initial MC event ranges built"); + unsigned int nROFRec = rofRecVec.size(); + + // Map hits for MC events + LOGP(info, "Map hits for MC events"); + mapHitsForMCEvents(hitTree, hitVecPool, mc2hitVec, mcEvMin, mcEvMax, nROFRec); + LOGP(info, "Mapped hits for MC events"); + + // Run cluster particle matching + auto pattIt = patternsPtr->cbegin(); + for (unsigned int iClus = 0; iClus < clusArr->size(); ++iClus) { + auto lab = (clusLabArr->getLabels(iClus))[0]; + const auto& c = (*clusArr)[iClus]; + // Ensure pattIt is advanced even if cluster is skipped + if (!lab.isValid() || lab.getSourceID() != 0 || !lab.isCorrect()) { + if (c.getPatternID() == CompCluster::InvalidPatternID) { + o2::itsmft::ClusterPattern::skipPattern(pattIt); + } + continue; + } + + int trackID{0}, evID{0}, srcID{0}; + bool fake{false}; + lab.get(trackID, evID, srcID, fake); + if (evID < 0 || static_cast(evID) >= info.size() || trackID < 0 || static_cast(trackID) >= info[evID].size()) { + if (c.getPatternID() == CompCluster::InvalidPatternID) { + o2::itsmft::ClusterPattern::skipPattern(pattIt); + } + continue; + } + UShort_t chipID = c.getSensorID(); + int layer = gman->getLayer(chipID); + bool isIB = layer < 3; + info[evID][trackID].clusters |= 1 << layer; + + o2::math_utils::Point3D clusterPos; + int clusterSize; + auto pattID = c.getPatternID(); + o2::itsmft::ClusterPattern patt; + if (!resolvePattern(c, pattIt, dict, isIB, patt)) { + continue; + } + clusterSize = patt.getNPixels(); + clusterPos = dict.getClusterCoordinates(c, patt, false); + + if (isIB) { + CurvedLocalToFlat(clusterPos, SegmentationIB(layer)); + } + + info[evID][trackID].clusterSize[layer] = clusterSize; + info[evID][trackID].clusterPattern[layer] = pattID; + info[evID][trackID].clusterLocX[layer] = clusterPos.X(); + info[evID][trackID].clusterLocZ[layer] = clusterPos.Z(); + + const auto& mc2hit = mc2hitVec[lab.getEventID()]; + const auto* hitArray = hitVecPool[lab.getEventID()]; + uint64_t key = (uint64_t(trackID) << 32) + c.getSensorID(); + auto hitIt = mc2hit.find(key); + if (hitIt == mc2hit.end()) + continue; + const auto& hit = (*hitArray)[hitIt->second]; + + auto hitLocSta = gman->getMatrixL2G(chipID) ^ hit.GetPosStart(); + auto hitLocEnd = gman->getMatrixL2G(chipID) ^ hit.GetPos(); + + if (isIB) { + CurvedLocalToFlat(hitLocSta, SegmentationIB(layer)); + CurvedLocalToFlat(hitLocEnd, SegmentationIB(layer)); + info[evID][trackID].hitLocX[layer] = 0.5f * (hitLocSta.X() + hitLocEnd.X()); + info[evID][trackID].hitLocY[layer] = 0.5f * (hitLocSta.Y() + hitLocEnd.Y()); + info[evID][trackID].hitLocZ[layer] = 0.5f * (hitLocSta.Z() + hitLocEnd.Z()); + } else { + auto x0 = hitLocSta.X(), dx = hitLocEnd.X() - x0; + auto y0 = hitLocSta.Y(), dy = hitLocEnd.Y() - y0; + auto z0 = hitLocSta.Z(), dz = hitLocEnd.Z() - z0; + auto r = (0.5f * (SegmentationOB::SensorLayerThickness - SegmentationOB::SensorLayerThicknessEff) - y0) / dy; + info[evID][trackID].hitLocX[layer] = x0 + r * dx; + info[evID][trackID].hitLocY[layer] = y0 + r * dy; + info[evID][trackID].hitLocZ[layer] = z0 + r * dz; + } + } + + LOGP(info, "Finished cluster-to-particle matching"); + + // The following part generates statistical histograms and outputs a TTree + int nb = 100; + double xbins[nb + 1], ptcutl = 0.01, ptcuth = 10.; + double a = std::log(ptcuth / ptcutl) / nb; + for (int i = 0; i <= nb; ++i) { + xbins[i] = ptcutl * std::exp(i * a); + } + auto* h_pt_num = new TH1D("h_pt_num", ";#it{p}_{T} (GeV/#it{c});Number of tracks", nb, xbins); + auto* h_pt_den = new TH1D("h_pt_den", ";#it{p}_{T} (GeV/#it{c});Number of generated primary particles", nb, xbins); + auto* h_pt_eff = new TEfficiency("h_pt_eff", "Tracking Efficiency;#it{p}_{T} (GeV/#it{c});Eff.", nb, xbins); + + auto* h_eta_num = new TH1D("h_eta_num", ";#it{#eta};Number of tracks", 60, -3, 3); + auto* h_eta_den = new TH1D("h_eta_den", ";#it{#eta};Number of generated particles", 60, -3, 3); + auto* h_eta_eff = new TEfficiency("h_eta_eff", "Tracking Efficiency;#it{#eta};Eff.", 60, -3, 3); + + auto* h_phi_num = new TH1D("h_phi_num", ";#varphi;Number of tracks", 360, 0., 2 * TMath::Pi()); + auto* h_phi_den = new TH1D("h_phi_den", ";#varphi;Number of generated particles", 360, 0., 2 * TMath::Pi()); + auto* h_phi_eff = new TEfficiency("h_phi_eff", "Tracking Efficiency;#varphi;Eff.", 360, 0., 2 * TMath::Pi()); + + auto* h_pt_fake = new TH1D("h_pt_fake", ";#it{p}_{T} (GeV/#it{c});Number of fake tracks", nb, xbins); + auto* h_pt_multifake = new TH1D("h_pt_multifake", ";#it{p}_{T} (GeV/#it{c});Number of multifake tracks", nb, xbins); + auto* h_pt_clones = new TH1D("h_pt_clones", ";#it{p}_{T} (GeV/#it{c});Number of cloned tracks", nb, xbins); + auto* h_dcaxy_vs_pt = new TH2D("h_dcaxy_vs_pt", ";#it{p}_{T} (GeV/#it{c});DCA_{xy} (#mum)", nb, xbins, 2000, -500., 500.); + auto* h_dcaxy_vs_eta = new TH2D("h_dcaxy_vs_eta", ";#it{#eta};DCA_{xy} (#mum)", 60, -3, 3, 2000, -500., 500.); + auto* h_dcaxy_vs_phi = new TH2D("h_dcaxy_vs_phi", ";#varphi;DCA_{xy} (#mum)", 360, 0., 2 * TMath::Pi(), 2000, -500., 500.); + auto* h_dcaz_vs_pt = new TH2D("h_dcaz_vs_pt", ";#it{p}_{T} (GeV/#it{c});DCA_{z} (#mum)", nb, xbins, 2000, -500., 500.); + auto* h_dcaz_vs_eta = new TH2D("h_dcaz_vs_eta", ";#it{#eta};DCA_{z} (#mum)", 60, -3, 3, 2000, -500., 500.); + auto* h_dcaz_vs_phi = new TH2D("h_dcaz_vs_phi", ";#varphi;DCA_{z} (#mum)", 360, 0., 2 * TMath::Pi(), 2000, -500., 500.); + auto* h_chi2 = new TH2D("h_chi2", ";#it{p}_{T} (GeV/#it{c});#chi^{2};Number of tracks", nb, xbins, 200, 0., 100.); + + for (auto& evInfo : info) { + for (auto& part : evInfo) { + if ((part.clusters & 0x7f) != 0x7f) { + // part.clusters != 0x3f && part.clusters != 0x3f << 1 && + // part.clusters != 0x1f && part.clusters != 0x1f << 1 && part.clusters + // != 0x1f << 2 && part.clusters != 0x0f && part.clusters != 0x0f << 1 + // && part.clusters != 0x0f << 2 && part.clusters != 0x0f << 3) { + continue; + } + if (!part.isPrimary) { + continue; + } + + h_pt_den->Fill(part.pt); + h_eta_den->Fill(part.eta); + h_phi_den->Fill(part.phi); + + if (part.isReco != 0u) { + h_pt_num->Fill(part.pt); + h_eta_num->Fill(part.eta); + h_phi_num->Fill(part.phi); + if (std::abs(part.eta) < 0.5) { + h_dcaxy_vs_pt->Fill(part.pt, part.dcaxy * 10000); + h_dcaz_vs_pt->Fill(part.pt, part.dcaz * 10000); + } + h_dcaz_vs_eta->Fill(part.eta, part.dcaz * 10000); + h_dcaxy_vs_eta->Fill(part.eta, part.dcaxy * 10000); + h_dcaxy_vs_phi->Fill(part.phi, part.dcaxy * 10000); + h_dcaz_vs_phi->Fill(part.phi, part.dcaz * 10000); + + h_chi2->Fill(part.pt, part.track.getChi2()); + + if (part.isReco > 1) { + for (int _i{0}; _i < part.isReco - 1; ++_i) { + h_pt_clones->Fill(part.pt); + } + } + } + if (part.isFake != 0u) { + h_pt_fake->Fill(part.pt); + if (part.isFake > 1) { + for (int _i{0}; _i < part.isFake - 1; ++_i) { + h_pt_multifake->Fill(part.pt); + } + } + } + } + } + + LOGP(info, "Streaming output TTree to file"); + TFile file("CorrTracksClusters.root", "recreate"); + TTree tree("ParticleInfo", "ParticleInfo"); + ParticleInfo pInfo; + tree.Branch("particle", &pInfo); + for (auto& event : info) { + for (auto& part : event) { + int nCl{0}; + for (unsigned int bit{0}; bit < sizeof(pInfo.clusters) * 8; ++bit) { + nCl += bool(part.clusters & (1 << bit)); + } + if (nCl < 3) { + continue; + } + pInfo = part; + tree.Fill(); + } + } + tree.Write(); + h_pt_num->Write(); + h_eta_num->Write(); + h_phi_num->Write(); + h_pt_den->Write(); + h_eta_den->Write(); + h_phi_den->Write(); + h_pt_multifake->Write(); + h_pt_fake->Write(); + h_dcaxy_vs_pt->Write(); + h_dcaz_vs_pt->Write(); + h_dcaxy_vs_eta->Write(); + h_dcaxy_vs_phi->Write(); + h_dcaz_vs_eta->Write(); + h_dcaz_vs_phi->Write(); + h_pt_clones->Write(); + h_chi2->Write(); + + h_pt_eff->SetTotalHistogram(*h_pt_den, ""); + h_pt_eff->SetPassedHistogram(*h_pt_num, ""); + h_pt_eff->SetTitle("Tracking Efficiency;#it{p}_{T} (GeV/#it{c});Eff."); + h_pt_eff->Write(); + + h_phi_eff->SetTotalHistogram(*h_phi_den, ""); + h_phi_eff->SetPassedHistogram(*h_phi_num, ""); + h_phi_eff->SetTitle("Tracking Efficiency;#varphi;Eff."); + h_phi_eff->Write(); + + h_eta_eff->SetTotalHistogram(*h_eta_den, ""); + h_eta_eff->SetPassedHistogram(*h_eta_num, ""); + h_eta_eff->SetTitle("Tracking Efficiency;#it{#eta};Eff."); + h_eta_eff->Write(); + + file.Close(); + LOGP(info, "Finished streaming output TTree to file"); + LOGP(info, "done."); +} diff --git a/Detectors/Upgrades/ITS3/macros/test/CreateDictionariesITS3.C b/Detectors/Upgrades/ITS3/macros/test/CreateDictionariesITS3.C index cc241afb3357a..76d7bf09de77f 100644 --- a/Detectors/Upgrades/ITS3/macros/test/CreateDictionariesITS3.C +++ b/Detectors/Upgrades/ITS3/macros/test/CreateDictionariesITS3.C @@ -52,7 +52,7 @@ #endif -void CreateDictionariesITS3(bool saveDeltas = false, +void CreateDictionariesITS3(bool saveDeltas = true, float probThreshold = 1e-6, std::string clusDictFile = "", std::string clusfile = "o2clus_its.root", @@ -94,7 +94,7 @@ void CreateDictionariesITS3(bool saveDeltas = false, TNtuple* nt = nullptr; if (saveDeltas) { fout = TFile::Open("CreateDictionaries.root", "recreate"); - nt = new TNtuple("nt", "hashes ntuple", "hash:dx:dz"); + nt = new TNtuple("nt", "hashes ntuple", "hash:layer:chipID:xhf:zhf:xcf:zcf:dx:dz:outlimDx:outlimDz"); } const o2::steer::DigitizationContext* digContext = nullptr; @@ -284,19 +284,25 @@ void CreateDictionariesITS3(bool saveDeltas = false, dZ = xyzLocM.Z() - locC.Z(); dX /= (ib) ? o2::its3::SegmentationMosaix::PitchRow : o2::itsmft::SegmentationAlpide::PitchRow; dZ /= (ib) ? o2::its3::SegmentationMosaix::PitchCol : o2::itsmft::SegmentationAlpide::PitchCol; - if (saveDeltas) { - nt->Fill(topology.getHash(), dX, dZ); - } + + float outLimitDx{-1}, outLimitDz{-1}; if (checkOutliers > 0.) { - if (bool bX = std::abs(dX) > topology.getRowSpan() * checkOutliers, bZ = std::abs(dZ) > topology.getColumnSpan() * checkOutliers; bX || bZ) { // ignore outlier + outLimitDx = topology.getRowSpan() * checkOutliers; + outLimitDz = topology.getColumnSpan() * checkOutliers; + bool isOutDx = std::abs(dX) > outLimitDx; + bool isOutDz = std::abs(dZ) > outLimitDz; + if (isOutDx || isOutDz) { // ignore outlier (ib) ? ++cOutliersIB : ++cOutliersOB; - LOGP(debug, "Ignored Value dX={} > {} * {} -> {}", dX, topology.getRowSpan(), checkOutliers, bX); - LOGP(debug, "Ignored Value dZ={} > {} * {} -> {}", dZ, topology.getColumnSpan(), checkOutliers, bZ); + LOGP(debug, "Ignored Value dX={} > {} * {} -> {}", dX, topology.getRowSpan(), checkOutliers, isOutDx); + LOGP(debug, "Ignored Value dZ={} > {} * {} -> {}", dZ, topology.getColumnSpan(), checkOutliers, isOutDz); dX = dZ = BuildTopologyDictionary::IgnoreVal; } else { (ib) ? ++cOkIB : ++cOkOB; } } + if (saveDeltas) { + nt->Fill(topology.getHash(), layer, chipID, xyzLocM.X(), xyzLocM.Z(), locC.X(), locC.Z(), dX, dZ, outLimitDx, outLimitDz); + } } } else { /* LOGP(info, " Failed to find MC hit entry for Tr: {} chipID: {}", trID, chipID); */ diff --git a/Detectors/Upgrades/ITS3/simulation/CMakeLists.txt b/Detectors/Upgrades/ITS3/simulation/CMakeLists.txt index 2fad72a96426d..8c4722012224d 100644 --- a/Detectors/Upgrades/ITS3/simulation/CMakeLists.txt +++ b/Detectors/Upgrades/ITS3/simulation/CMakeLists.txt @@ -15,8 +15,11 @@ o2_add_library(ITS3Simulation src/DescriptorInnerBarrelITS3.cxx src/Digitizer.cxx src/DigiParams.cxx + src/ITS3DPLDigitizerParam.cxx + src/ChipDigitsContainer.cxx + src/ChipSimResponse.cxx PUBLIC_LINK_LIBRARIES O2::SimulationDataFormat - O2::ITSBase O2::ITSMFTSimulation + O2::ITSBase O2::ITSMFTSimulation O2::ITSMFTBase ROOT::Physics) o2_target_root_dictionary(ITS3Simulation @@ -25,6 +28,9 @@ o2_target_root_dictionary(ITS3Simulation include/ITS3Simulation/DescriptorInnerBarrelITS3.h include/ITS3Simulation/Digitizer.h include/ITS3Simulation/DigiParams.h + include/ITS3Simulation/ITS3DPLDigitizerParam.h + include/ITS3Simulation/ChipDigitsContainer.h + include/ITS3Simulation/ChipSimResponse.h ) o2_data_file(COPY data DESTINATION Detectors/ITS3/simulation) diff --git a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ChipDigitsContainer.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ChipDigitsContainer.h new file mode 100644 index 0000000000000..0c9627fe412c3 --- /dev/null +++ b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ChipDigitsContainer.h @@ -0,0 +1,59 @@ +// 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_CHIPDIGITSCONTAINER_ +#define ALICEO2_ITS3_CHIPDIGITSCONTAINER_ + +#include "ITSMFTBase/SegmentationAlpide.h" // Base class in o2::itsmft namespace +#include "ITSMFTSimulation/ChipDigitsContainer.h" // Base class in o2::itsmft namespace +#include "ITS3Base/SegmentationMosaix.h" // OB segmentation implementation +#include "ITS3Base/SpecsV2.h" // Provides SpecsV2::isDetITS3() interface +#include "ITS3Simulation/DigiParams.h" // ITS3-specific DigiParams interface +#include + +namespace o2::its3 +{ + +class ChipDigitsContainer : public o2::itsmft::ChipDigitsContainer +{ + private: + bool innerBarrel; ///< true if the chip belongs to the inner barrel (IB), false if outer barrel (OB) + int maxRows; ///< maximum number of rows + int maxCols; ///< maximum number of columns + + public: + explicit ChipDigitsContainer(UShort_t idx = 0); + + using SegmentationIB = SegmentationMosaix; + using SegmentationOB = o2::itsmft::SegmentationAlpide; + + /// Returns whether the chip is in the inner barrel (IB) + void setChipIndex(UShort_t idx) + { + o2::itsmft::ChipDigitsContainer::setChipIndex(idx); + innerBarrel = constants::detID::isDetITS3(getChipIndex()); + maxRows = innerBarrel ? SegmentationIB::NRows : SegmentationOB::NRows; + maxCols = innerBarrel ? SegmentationIB::NCols : SegmentationOB::NCols; + } + + int getMaxRows() const { return maxRows; } + int getMaxCols() const { return maxCols; } + bool isIB() const; + /// Adds noise digits, deleted the one using the itsmft::DigiParams interface + void addNoise(UInt_t rofMin, UInt_t rofMax, const o2::itsmft::DigiParams* params, int maxRows = o2::itsmft::SegmentationAlpide::NRows, int maxCols = o2::itsmft::SegmentationAlpide::NCols) = delete; + void addNoise(UInt_t rofMin, UInt_t rofMax, const o2::its3::DigiParams* params); + + ClassDefNV(ChipDigitsContainer, 1); +}; + +} // namespace o2::its3 + +#endif // ALICEO2_ITS3_CHIPDIGITSCONTAINER_ \ No newline at end of file diff --git a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ChipSimResponse.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ChipSimResponse.h new file mode 100644 index 0000000000000..f96fde9fb0d55 --- /dev/null +++ b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ChipSimResponse.h @@ -0,0 +1,41 @@ +// 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_ITS3SIMULATION_CHIPSIMRESPONSE_H +#define ALICEO2_ITS3SIMULATION_CHIPSIMRESPONSE_H + +#include "ITSMFTSimulation/AlpideSimResponse.h" + +namespace o2 +{ +namespace its3 +{ + +class ChipSimResponse : public o2::itsmft::AlpideSimResponse +{ + public: + ChipSimResponse() = default; + ChipSimResponse(const ChipSimResponse& other) = default; + + float getRespCentreDep() const { return mRespCentreDep; } + void computeCentreFromData(); + void initData(int tableNumber, std::string dataPath, const bool quiet = true); + + private: + float mRespCentreDep = 0.f; + + ClassDef(ChipSimResponse, 1); +}; + +} // namespace its3 +} // namespace o2 + +#endif // ALICEO2_ITS3SIMULATION_CHIPSIMRESPONSE_H \ No newline at end of file diff --git a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/DigiParams.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/DigiParams.h index eca0a71949ba7..5764dfbd7d593 100644 --- a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/DigiParams.h +++ b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/DigiParams.h @@ -13,21 +13,43 @@ #define ITS3_DIGIPARAMS_H #include "ITSMFTSimulation/DigiParams.h" +#include "ITS3Simulation/ChipSimResponse.h" namespace o2::its3 { class DigiParams final : public o2::itsmft::DigiParams { + private: + float mIBNoisePerPixel = 1.e-8; + int mIBChargeThreshold = 150; ///< charge threshold in Nelectrons + int mIBMinChargeToAccount = 15; ///< minimum charge contribution to account + int mIBNSimSteps = 18; ///< number of steps in response simulation + float mIBNSimStepsInv = 0; ///< its inverse + public: + DigiParams(); + + void setIBNoisePerPixel(float v) { mIBNoisePerPixel = v; } + float getIBNoisePerPixel() const { return mIBNoisePerPixel; } + + void setIBChargeThreshold(int v, float frac2Account = 0.1); + int getIBChargeThreshold() const { return mIBChargeThreshold; } + + void setIBNSimSteps(int v); + int getIBNSimSteps() const { return mIBNSimSteps; } + float getIBNSimStepsInv() const { return mIBNSimStepsInv; } + + int getIBMinChargeToAccount() const { return mIBMinChargeToAccount; } + const o2::itsmft::AlpideSimResponse* getAlpSimResponse() const = delete; void setAlpSimResponse(const o2::itsmft::AlpideSimResponse* par) = delete; const o2::itsmft::AlpideSimResponse* getOBSimResponse() const { return mOBSimResponse; } void setOBSimResponse(const o2::itsmft::AlpideSimResponse* response) { mOBSimResponse = response; } - const o2::itsmft::AlpideSimResponse* getIBSimResponse() const { return mIBSimResponse; } - void setIBSimResponse(const o2::itsmft::AlpideSimResponse* response) { mIBSimResponse = response; } + o2::its3::ChipSimResponse* getIBSimResponse() const { return mIBSimResponse; } + void setIBSimResponse(o2::its3::ChipSimResponse* response); bool hasResponseFunctions() const { return mIBSimResponse != nullptr && mOBSimResponse != nullptr; } @@ -35,7 +57,7 @@ class DigiParams final : public o2::itsmft::DigiParams private: const o2::itsmft::AlpideSimResponse* mOBSimResponse = nullptr; //!< pointer to external response - const o2::itsmft::AlpideSimResponse* mIBSimResponse = nullptr; //!< pointer to external response + o2::its3::ChipSimResponse* mIBSimResponse = nullptr; //!< pointer to external response ClassDef(DigiParams, 1); }; diff --git a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h index 8d0f06a27343b..edc5583c03d5a 100644 --- a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h +++ b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h @@ -21,12 +21,13 @@ #include "Rtypes.h" #include "TObject.h" -#include "ITSMFTSimulation/ChipDigitsContainer.h" #include "ITSMFTSimulation/AlpideSimResponse.h" #include "ITSMFTSimulation/Hit.h" #include "ITSBase/GeometryTGeo.h" #include "ITS3Base/SegmentationMosaix.h" #include "ITS3Simulation/DigiParams.h" +#include "ITS3Simulation/ChipDigitsContainer.h" +#include "ITS3Simulation/ChipSimResponse.h" #include "DataFormatsITSMFT/Digit.h" #include "DataFormatsITSMFT/ROFRecord.h" #include "CommonDataFormat/InteractionRecord.h" @@ -78,7 +79,7 @@ class Digitizer : public TObject private: void processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID, int srcID); - void registerDigits(o2::itsmft::ChipDigitsContainer& chip, uint32_t roFrame, float tInROF, int nROF, + void registerDigits(o2::its3::ChipDigitsContainer& chip, uint32_t roFrame, float tInROF, int nROF, uint16_t row, uint16_t col, int nEle, o2::MCCompLabel& lbl); ExtraDig* getExtraDigBuffer(uint32_t roFrame) @@ -108,7 +109,7 @@ class Digitizer : public TObject static constexpr std::array mIBSegmentations{0, 1, 2}; - o2::itsmft::AlpideSimResponse* mSimRespIB = nullptr; // simulated response for IB + o2::its3::ChipSimResponse* mSimRespIB = nullptr; // simulated response for IB o2::itsmft::AlpideSimResponse* mSimRespOB = nullptr; // simulated response for OB bool mSimRespIBOrientation{false}; // wether the orientation in the IB response function is flipped float mSimRespIBShift{0.f}; // adjusting the Y-shift in the IB response function to match sensor local coord. @@ -118,7 +119,7 @@ class Digitizer : public TObject const o2::its::GeometryTGeo* mGeometry = nullptr; ///< ITS3 geometry - std::vector mChips; ///< Array of chips digits containers + std::vector mChips; ///< Array of chips digits containers std::deque> mExtraBuff; ///< burrer (per roFrame) for extra digits std::vector* mDigits = nullptr; //! output digits diff --git a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ITS3DPLDigitizerParam.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ITS3DPLDigitizerParam.h new file mode 100644 index 0000000000000..3192f73fb8f79 --- /dev/null +++ b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ITS3DPLDigitizerParam.h @@ -0,0 +1,32 @@ +// 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_ITS3DPLDIGITIZERPARAM_H_ +#define ALICEO2_ITS3DPLDIGITIZERPARAM_H_ + +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/ConfigurableParamHelper.h" + +namespace o2::its3 +{ + +struct ITS3DPLDigitizerParam : public o2::conf::ConfigurableParamHelper { + float IBNoisePerPixel = 1.e-8; ///< MOSAIX Noise per channel + int IBChargeThreshold = 150; ///< charge threshold in Nelectrons for IB + int IBMinChargeToAccount = 15; ///< minimum charge contribution to account for IB + int nIBSimSteps = 18; ///< number of steps in response for IB + + O2ParamDef(ITS3DPLDigitizerParam, "ITS3DPLDigitizerParam"); +}; + +} // namespace o2::its3 + +#endif \ No newline at end of file diff --git a/Detectors/Upgrades/ITS3/simulation/src/ChipDigitsContainer.cxx b/Detectors/Upgrades/ITS3/simulation/src/ChipDigitsContainer.cxx new file mode 100644 index 0000000000000..0611f7002f160 --- /dev/null +++ b/Detectors/Upgrades/ITS3/simulation/src/ChipDigitsContainer.cxx @@ -0,0 +1,63 @@ +// 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 "ITS3Simulation/ChipDigitsContainer.h" + +namespace o2 +{ +namespace its3 +{ + +ChipDigitsContainer::ChipDigitsContainer(UShort_t idx) + : o2::itsmft::ChipDigitsContainer(idx) {} + +bool ChipDigitsContainer::isIB() const +{ + return innerBarrel; +} + +void ChipDigitsContainer::addNoise(UInt_t rofMin, UInt_t rofMax, const o2::its3::DigiParams* params) +{ + UInt_t row = 0; + UInt_t col = 0; + Int_t nhits = 0; + constexpr float ns2sec = 1e-9; + float mean = 0.f; + int nel = 0; + + if (isIB()) { + // Inner barrel: use ITS3-specific noise interface with OB segmentation. + mean = params->getIBNoisePerPixel() * SegmentationOB::NPixels; + nel = static_cast(params->getIBChargeThreshold() * 1.1); + } else { + // Outer barrel: use base class noise interface with IB segmentation. + mean = params->getNoisePerPixel() * SegmentationIB::NPixels; + nel = static_cast(params->getChargeThreshold() * 1.1); + } + + for (UInt_t rof = rofMin; rof <= rofMax; ++rof) { + nhits = gRandom->Poisson(mean); + for (Int_t i = 0; i < nhits; ++i) { + row = gRandom->Integer(maxRows); + col = gRandom->Integer(maxCols); + if (mNoiseMap && mNoiseMap->isNoisy(mChipIndex, row, col)) + continue; + if (mDeadChanMap && mDeadChanMap->isNoisy(mChipIndex, row, col)) + continue; + auto key = getOrderingKey(rof, row, col); + if (!findDigit(key)) + addDigit(key, rof, row, col, nel, o2::MCCompLabel(true)); + } + } +} + +} // namespace its3 +} // namespace o2 \ No newline at end of file diff --git a/Detectors/Upgrades/ITS3/simulation/src/ChipSimResponse.cxx b/Detectors/Upgrades/ITS3/simulation/src/ChipSimResponse.cxx new file mode 100644 index 0000000000000..1c482983f0d0a --- /dev/null +++ b/Detectors/Upgrades/ITS3/simulation/src/ChipSimResponse.cxx @@ -0,0 +1,62 @@ +// 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 "ITS3Simulation/ChipSimResponse.h" +#include +#include + +using namespace o2::its3; + +ClassImp(o2::its3::ChipSimResponse); + +void ChipSimResponse::initData(int tableNumber, std::string dataPath, const bool quiet) +{ + AlpideSimResponse::initData(tableNumber, dataPath, quiet); + computeCentreFromData(); +} + +void ChipSimResponse::computeCentreFromData() +{ + std::vector zVec, qVec; + const int npix = o2::itsmft::AlpideRespSimMat::getNPix(); + + for (int iz = 0; iz < mNBinDpt; ++iz) { + size_t bin = iz + mNBinDpt * (0 + mNBinRow * 0); + const auto& mat = mData[bin]; + float val = mat.getValue(npix / 2, npix / 2); + float gz = mDptMin + iz / mStepInvDpt; + zVec.push_back(gz); + qVec.push_back(val); + } + + std::vector> zqPairs; + for (size_t i = 0; i < zVec.size(); ++i) { + zqPairs.emplace_back(zVec[i], qVec[i]); + } + std::sort(zqPairs.begin(), zqPairs.end()); + zVec.clear(); + qVec.clear(); + for (auto& p : zqPairs) { + zVec.push_back(p.first); + qVec.push_back(p.second); + } + + float intQ = 0.f, intZQ = 0.f; + for (size_t i = 0; i + 1 < zVec.size(); ++i) { + float z0 = zVec[i], z1 = zVec[i + 1]; + float q0 = qVec[i], q1 = qVec[i + 1]; + float dz = z1 - z0; + intQ += 0.5f * (q0 + q1) * dz; + intZQ += 0.5f * (z0 * q0 + z1 * q1) * dz; + } + + mRespCentreDep = (intQ > 0.f) ? intZQ / intQ : 0.f; +} diff --git a/Detectors/Upgrades/ITS3/simulation/src/DigiParams.cxx b/Detectors/Upgrades/ITS3/simulation/src/DigiParams.cxx index a9f17a544b3c4..afa02ec44741d 100644 --- a/Detectors/Upgrades/ITS3/simulation/src/DigiParams.cxx +++ b/Detectors/Upgrades/ITS3/simulation/src/DigiParams.cxx @@ -14,27 +14,67 @@ #include "Framework/Logger.h" #include "ITS3Simulation/DigiParams.h" +#include ClassImp(o2::its3::DigiParams); namespace o2::its3 { +DigiParams::DigiParams() +{ + // make sure the defaults are consistent + setIBNSimSteps(mIBNSimSteps); +} + +void DigiParams::setIBNSimSteps(int v) +{ + // set number of sampling steps in silicon + mIBNSimSteps = v > 0 ? v : 1; + mIBNSimStepsInv = 1.f / mIBNSimSteps; +} + +void DigiParams::setIBChargeThreshold(int v, float frac2Account) +{ + // set charge threshold for digits creation and its fraction to account + // contribution from single hit + mIBChargeThreshold = v; + mIBMinChargeToAccount = v * frac2Account; + if (mIBMinChargeToAccount < 0 || mIBMinChargeToAccount > mIBChargeThreshold) { + mIBMinChargeToAccount = mIBChargeThreshold; + } + LOG(info) << "Set Mosaix charge threshold to " << mIBChargeThreshold + << ", single hit will be accounted from " << mIBMinChargeToAccount + << " electrons"; +} + void DigiParams::print() const { // print settings - LOGF(info, "ITS3 DigiParams settings:"); - LOGF(info, "Continuous readout : %s", isContinuous() ? "ON" : "OFF"); - LOGF(info, "Readout Frame Length(ns) : %f", getROFrameLength()); - LOGF(info, "Strobe delay (ns) : %f", getStrobeDelay()); - LOGF(info, "Strobe length (ns) : %f", getStrobeLength()); - LOGF(info, "Threshold (N electrons) : %d", getChargeThreshold()); - LOGF(info, "Min N electrons to account : %d", getMinChargeToAccount()); - LOGF(info, "Number of charge sharing steps : %d", getNSimSteps()); - LOGF(info, "ELoss to N electrons factor : %e", getEnergyToNElectrons()); - LOGF(info, "Noise level per pixel : %e", getNoisePerPixel()); - LOGF(info, "Charge time-response:\n"); + printf("ITS3 DigiParams settings:\n"); + printf("Continuous readout : %s\n", isContinuous() ? "ON" : "OFF"); + printf("Readout Frame Length(ns) : %f\n", getROFrameLength()); + printf("Strobe delay (ns) : %f\n", getStrobeDelay()); + printf("Strobe length (ns) : %f\n", getStrobeLength()); + printf("IB Threshold (N electrons) : %d\n", getIBChargeThreshold()); + printf("OB Threshold (N electrons) : %d\n", getChargeThreshold()); + printf("Min N electrons to account for IB : %d\n", getIBMinChargeToAccount()); + printf("Min N electrons to account for OB : %d\n", getMinChargeToAccount()); + printf("Number of charge sharing steps of IB : %d\n", getIBNSimSteps()); + printf("Number of charge sharing steps of OB : %d\n", getNSimSteps()); + printf("ELoss to N electrons factor : %e\n", getEnergyToNElectrons()); + printf("Noise level per pixel of IB : %e\n", getIBNoisePerPixel()); + printf("Noise level per pixel of OB : %e\n", getNoisePerPixel()); + printf("Charge time-response:\n"); getSignalShape().print(); } +void DigiParams::setIBSimResponse(o2::its3::ChipSimResponse* response) +{ + mIBSimResponse = response; + if (mIBSimResponse) { + mIBSimResponse->computeCentreFromData(); + } +} + } // namespace o2::its3 diff --git a/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx index 3c75bf3e8f680..1d1d15a91f89b 100644 --- a/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx @@ -27,7 +27,8 @@ #include using o2::itsmft::Hit; -using SegmentationAlpide = o2::itsmft::SegmentationAlpide; +using SegmentationOB = o2::itsmft::SegmentationAlpide; +using SegmentationIB = o2::its3::SegmentationMosaix; using o2::itsmft::AlpideRespSimMat; using o2::itsmft::PreDigit; @@ -46,8 +47,8 @@ void Digitizer::init() } if (!mParams.hasResponseFunctions()) { - auto loadSetResponseFunc = [&](const char* name, const char* fileIB, const char* nameIB, const char* fileOB, const char* nameOB) { - LOGP(info, "Loading response function for {}: IB={}:{} ; OB={}:{}", name, nameIB, fileIB, nameOB, fileOB); + auto loadSetResponseFunc = [&](const char* fileIB, const char* nameIB, const char* fileOB, const char* nameOB) { + LOGP(info, "Loading response function IB={}:{} ; OB={}:{}", nameIB, fileIB, nameOB, fileOB); auto fIB = TFile::Open(fileIB, "READ"); if (!fIB || fIB->IsZombie() || !fIB->IsOpen()) { LOGP(fatal, "Cannot open file {}", fileIB); @@ -56,7 +57,7 @@ void Digitizer::init() if (!fOB || fOB->IsZombie() || !fOB->IsOpen()) { LOGP(fatal, "Cannot open file {}", fileOB); } - mParams.setIBSimResponse(mSimRespIB = fIB->Get(nameIB)); + mParams.setIBSimResponse(mSimRespIB = fIB->Get(nameIB)); mParams.setOBSimResponse(mSimRespOB = fOB->Get(nameOB)); fIB->Close(); fOB->Close(); @@ -64,25 +65,27 @@ void Digitizer::init() if (const auto& func = ITS3Params::Instance().chipResponseFunction; func == "Alpide") { constexpr const char* responseFile = "$(O2_ROOT)/share/Detectors/ITSMFT/data/AlpideResponseData/AlpideResponseData.root"; - loadSetResponseFunc("Alpide", responseFile, "response0", responseFile, "response1"); - mSimRespIBShift = mSimRespIB->getDepthMax() - SegmentationMosaix::SensorLayerThickness / 2.f + 10.e-4f; - mSimRespOBShift = mSimRespOB->getDepthMax() - SegmentationAlpide::SensorLayerThickness / 2.f; + loadSetResponseFunc(responseFile, "response0", responseFile, "response0"); + mSimRespIBScaleX = o2::itsmft::SegmentationAlpide::PitchRow / SegmentationIB::PitchRow; + mSimRespIBScaleZ = o2::itsmft::SegmentationAlpide::PitchCol / SegmentationIB::PitchCol; } else if (func == "APTS") { constexpr const char* responseFileIB = "$(O2_ROOT)/share/Detectors/Upgrades/ITS3/data/ITS3ChipResponseData/APTSResponseData.root"; constexpr const char* responseFileOB = "$(O2_ROOT)/share/Detectors/ITSMFT/data/AlpideResponseData/AlpideResponseData.root"; - loadSetResponseFunc("APTS", responseFileIB, "response1", responseFileOB, "response1"); - mSimRespIBShift = mSimRespIB->getDepthMax() + (float)constants::pixelarray::pixels::apts::responseYShift; - mSimRespOBShift = mSimRespOB->getDepthMax() - SegmentationAlpide::SensorLayerThickness / 2.f; - mSimRespIBScaleX = 0.5f * constants::pixelarray::pixels::apts::pitchX / SegmentationMosaix::PitchRow; - mSimRespIBScaleZ = 0.5f * constants::pixelarray::pixels::apts::pitchZ / SegmentationMosaix::PitchCol; + loadSetResponseFunc(responseFileIB, "response1", responseFileOB, "response0"); + mSimRespIBScaleX = constants::pixelarray::pixels::apts::pitchX / SegmentationIB::PitchRow; + mSimRespIBScaleZ = constants::pixelarray::pixels::apts::pitchZ / SegmentationIB::PitchCol; mSimRespIBOrientation = true; } else { LOGP(fatal, "ResponseFunction '{}' not implemented!", func); } + mSimRespIBShift = mSimRespIB->getDepthMax() - constants::silicon::thickness / 2.f; + mSimRespOBShift = mSimRespOB->getDepthMax() - SegmentationOB::SensorLayerThickness / 2.f; } + mParams.print(); - LOGP(info, "IBShift = {} ; OBShift = {}", mSimRespIBShift, mSimRespOBShift); - LOGP(info, "IB-Scale: X={} ; Z={}", mSimRespIBScaleX, mSimRespIBScaleZ); + LOGP(info, "IB shift = {} ; OB shift = {}", mSimRespIBShift, mSimRespOBShift); + LOGP(info, "IB pixel scale on x = {} ; z = {}", mSimRespIBScaleX, mSimRespIBScaleZ); + LOGP(info, "IB response orientation: {}", mSimRespIBOrientation ? "flipped" : "normal"); mIRFirstSampledTF = o2::raw::HBFUtils::Instance().getFirstSampledTFIR(); } @@ -173,11 +176,7 @@ void Digitizer::fillOutputContainer(uint32_t frameLast) auto& extra = *(mExtraBuff.front().get()); for (size_t iChip{0}; iChip < mChips.size(); ++iChip) { auto& chip = mChips[iChip]; - if (constants::detID::isDetITS3(iChip)) { // Check if this is a chip of ITS3 - chip.addNoise(mROFrameMin, mROFrameMin, &mParams, SegmentationMosaix::NRows, SegmentationMosaix::NCols); - } else { - chip.addNoise(mROFrameMin, mROFrameMin, &mParams); - } + chip.addNoise(mROFrameMin, mROFrameMin, &mParams); auto& buffer = chip.getPreDigits(); if (buffer.empty()) { continue; @@ -190,7 +189,7 @@ void Digitizer::fillOutputContainer(uint32_t frameLast) break; // is the digit ROFrame from the key > the max requested frame } auto& preDig = iter->second; // preDigit - if (preDig.charge >= mParams.getChargeThreshold()) { + if (preDig.charge >= (chip.isIB() ? mParams.getIBChargeThreshold() : mParams.getChargeThreshold())) { int digID = mDigits->size(); mDigits->emplace_back(chip.getChipIndex(), preDig.row, preDig.col, preDig.charge); mMCLabels->addElement(digID, preDig.labelRef.label); @@ -257,16 +256,15 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID } // here we start stepping in the depth of the sensor to generate charge diffision - float nStepsInv = mParams.getNSimStepsInv(); - int nSteps = mParams.getNSimSteps(); int detID{hit.GetDetectorID()}; int layer = mGeometry->getLayer(detID); const auto& matrix = mGeometry->getMatrixL2G(detID); - bool innerBarrel{layer < 3}; + int nSteps = chip.isIB() ? mParams.getIBNSimSteps() : mParams.getNSimSteps(); + float nStepsInv = chip.isIB() ? mParams.getIBNSimStepsInv() : mParams.getNSimStepsInv(); math_utils::Vector3D xyzLocS, xyzLocE; xyzLocS = matrix ^ (hit.GetPosStart()); // Global hit coordinates to local detector coordinates xyzLocE = matrix ^ (hit.GetPos()); - if (innerBarrel) { + if (chip.isIB()) { // transform the point on the curved surface to a flat one float xFlatE{0.f}, yFlatE{0.f}, xFlatS{0.f}, yFlatS{0.f}; mIBSegmentations[layer].curvedToFlat(xyzLocS.X(), xyzLocS.Y(), xFlatS, yFlatS); @@ -284,7 +282,7 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID 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) { + if (chip.isIB()) { // get entrance pixel row and col while (!mIBSegmentations[layer].localToDetector(xyzLocS.X(), xyzLocS.Z(), rowS, colS)) { // guard-ring ? if (++nSkip >= nSteps) { @@ -301,14 +299,14 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID } } else { // get entrance pixel row and col - while (!SegmentationAlpide::localToDetector(xyzLocS.X(), xyzLocS.Z(), rowS, colS)) { // guard-ring ? + while (!SegmentationOB::localToDetector(xyzLocS.X(), xyzLocS.Z(), rowS, colS)) { // guard-ring ? if (++nSkip >= nSteps) { return; // did not enter to sensitive matrix } xyzLocS += step; } // get exit pixel row and col - while (!SegmentationAlpide::localToDetector(xyzLocE.X(), xyzLocE.Z(), rowE, colE)) { // guard-ring ? + while (!SegmentationOB::localToDetector(xyzLocE.X(), xyzLocE.Z(), rowE, colE)) { // guard-ring ? if (++nSkip >= nSteps) { return; // did not enter to sensitive matrix } @@ -327,8 +325,8 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID rowE += AlpideRespSimMat::NPix / 2; rowS = std::max(rowS, 0); - const int maxNrows{innerBarrel ? SegmentationMosaix::NRows : SegmentationAlpide::NRows}; - const int maxNcols{innerBarrel ? SegmentationMosaix::NCols : SegmentationAlpide::NCols}; + const int maxNrows{chip.isIB() ? SegmentationIB::NRows : SegmentationOB::NRows}; + const int maxNcols{chip.isIB() ? SegmentationIB::NCols : SegmentationOB::NCols}; rowE = std::min(rowE, maxNrows - 1); colS -= AlpideRespSimMat::NPix / 2; @@ -352,22 +350,22 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID // take into account that the AlpideSimResponse depth defintion has different min/max boundaries // although the max should coincide with the surface of the epitaxial layer, which in the chip // local coordinates has Y = +SensorLayerThickness/2 - xyzLocS.SetY(xyzLocS.Y() + ((innerBarrel) ? mSimRespIBShift : mSimRespOBShift)); + xyzLocS.SetY(xyzLocS.Y() + ((chip.isIB()) ? mSimRespIBShift : mSimRespOBShift)); // collect charge in evey pixel which might be affected by the hit for (int iStep = nSteps; iStep--;) { // Get the pixel ID - if (innerBarrel) { + if (chip.isIB()) { mIBSegmentations[layer].localToDetector(xyzLocS.X(), xyzLocS.Z(), row, col); } else { - SegmentationAlpide::localToDetector(xyzLocS.X(), xyzLocS.Z(), row, col); + SegmentationOB::localToDetector(xyzLocS.X(), xyzLocS.Z(), row, col); } if (row != rowPrev || col != colPrev) { // update pixel and coordinates of its center - if (innerBarrel) { + if (chip.isIB()) { if (!mIBSegmentations[layer].detectorToLocal(row, col, cRowPix, cColPix)) { continue; } - } else if (!SegmentationAlpide::detectorToLocal(row, col, cRowPix, cColPix)) { + } else if (!SegmentationOB::detectorToLocal(row, col, cRowPix, cColPix)) { continue; // should not happen } rowPrev = row; @@ -377,13 +375,13 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID // note that response needs coordinates along column row (locX) (locZ) then depth (locY) float rowMax{}, colMax{}; const AlpideRespSimMat* rspmat{nullptr}; - if (innerBarrel) { - rowMax = 0.5f * SegmentationMosaix::PitchRow; - colMax = 0.5f * SegmentationMosaix::PitchCol; + if (chip.isIB()) { + rowMax = 0.5f * SegmentationIB::PitchRow * mSimRespIBScaleX; + colMax = 0.5f * SegmentationIB::PitchCol * mSimRespIBScaleZ; rspmat = mSimRespIB->getResponse(mSimRespIBScaleX * (xyzLocS.X() - cRowPix), mSimRespIBScaleZ * (xyzLocS.Z() - cColPix), xyzLocS.Y(), flipRow, flipCol, rowMax, colMax); } else { - rowMax = 0.5f * SegmentationAlpide::PitchRow; - colMax = 0.5f * SegmentationAlpide::PitchCol; + rowMax = 0.5f * SegmentationOB::PitchRow; + colMax = 0.5f * SegmentationOB::PitchCol; rspmat = mSimRespOB->getResponse(xyzLocS.X() - cRowPix, xyzLocS.Z() - cColPix, xyzLocS.Y(), flipRow, flipCol, rowMax, colMax); } @@ -402,7 +400,7 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID if (colDest < 0 || colDest >= colSpan) { continue; } - respMatrix[rowDest][colDest] += rspmat->getValue(irow, icol, ((innerBarrel && mSimRespIBOrientation) ? !flipRow : flipRow), flipCol); + respMatrix[rowDest][colDest] += rspmat->getValue(irow, icol, ((chip.isIB() && mSimRespIBOrientation) ? !flipRow : flipRow), flipCol); } } } @@ -419,7 +417,7 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID } int nEle = gRandom->Poisson(nElectrons * nEleResp); // total charge in given pixel // ignore charge which have no chance to fire the pixel - if (nEle < mParams.getMinChargeToAccount()) { + if (nEle < (chip.isIB() ? mParams.getIBChargeThreshold() : mParams.getChargeThreshold())) { continue; } uint16_t colIS = icol + colS; @@ -428,7 +426,7 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID } } -void Digitizer::registerDigits(o2::itsmft::ChipDigitsContainer& chip, uint32_t roFrame, float tInROF, int nROF, +void Digitizer::registerDigits(o2::its3::ChipDigitsContainer& chip, uint32_t roFrame, float tInROF, int nROF, uint16_t row, uint16_t col, int nEle, o2::MCCompLabel& lbl) { // Register digits for given pixel, accounting for the possible signal contribution to @@ -442,7 +440,7 @@ void Digitizer::registerDigits(o2::itsmft::ChipDigitsContainer& chip, uint32_t r tStrobe += mParams.getROFrameLength(); // for the next ROF // discard too small contributions, they have no chance to produce a digit - if (nEleROF < mParams.getMinChargeToAccount()) { + if (nEleROF < (chip.isIB() ? mParams.getIBChargeThreshold() : mParams.getChargeThreshold())) { continue; } if (roFr > mEventROFrameMax) { diff --git a/Detectors/Upgrades/ITS3/simulation/src/ITS3DPLDigitizerParam.cxx b/Detectors/Upgrades/ITS3/simulation/src/ITS3DPLDigitizerParam.cxx new file mode 100644 index 0000000000000..69314b8a0be9b --- /dev/null +++ b/Detectors/Upgrades/ITS3/simulation/src/ITS3DPLDigitizerParam.cxx @@ -0,0 +1,14 @@ +// 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 "ITS3Simulation/ITS3DPLDigitizerParam.h" + +O2ParamImpl(o2::its3::ITS3DPLDigitizerParam) \ No newline at end of file diff --git a/Detectors/Upgrades/ITS3/simulation/src/ITS3SimulationLinkDef.h b/Detectors/Upgrades/ITS3/simulation/src/ITS3SimulationLinkDef.h index fca3f5d63c2c4..921512193f98b 100644 --- a/Detectors/Upgrades/ITS3/simulation/src/ITS3SimulationLinkDef.h +++ b/Detectors/Upgrades/ITS3/simulation/src/ITS3SimulationLinkDef.h @@ -20,5 +20,8 @@ #pragma link C++ class o2::its3::DescriptorInnerBarrelITS3 + ; #pragma link C++ class o2::its3::DigiParams + ; #pragma link C++ class o2::its3::Digitizer + ; +#pragma link C++ class o2::its3::ITS3DPLDigitizerParam + ; +#pragma link C++ class o2::its3::ChipDigitsContainer + ; +#pragma link C++ class o2::its3::ChipSimResponse + ; #endif diff --git a/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.cxx b/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.cxx index 27f876f7bc24b..af0af091d40e8 100644 --- a/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.cxx +++ b/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.cxx @@ -27,6 +27,7 @@ #include "DataFormatsITSMFT/ROFRecord.h" #include "ITS3Simulation/Digitizer.h" #include "ITSMFTSimulation/DPLDigitizerParam.h" +#include "ITS3Simulation/ITS3DPLDigitizerParam.h" #include "ITSMFTBase/DPLAlpideParam.h" #include "ITSBase/GeometryTGeo.h" #include "ITS3Base/ITS3Params.h" @@ -216,6 +217,7 @@ class ITS3DPLDigitizerTask : BaseDPLDigitizer mDigitizer.setGeometry(geom); const auto& dopt = o2::itsmft::DPLDigitizerParam::Instance(); + const auto& doptIB = o2::its3::ITS3DPLDigitizerParam::Instance(); pc.inputs().get*>("ITS_alppar"); const auto& aopt = o2::itsmft::DPLAlpideParam::Instance(); digipar.setContinuous(dopt.continuous); @@ -238,6 +240,11 @@ class ITS3DPLDigitizerTask : BaseDPLDigitizer digipar.setTimeOffset(dopt.timeOffset); digipar.setNSimSteps(dopt.nSimSteps); + // ITS3 inner barrel specific parameters + digipar.setIBChargeThreshold(doptIB.IBChargeThreshold); + digipar.setIBNSimSteps(doptIB.nIBSimSteps); + digipar.setIBNoisePerPixel(doptIB.IBNoisePerPixel); + mROMode = digipar.isContinuous() ? o2::parameters::GRPObject::CONTINUOUS : o2::parameters::GRPObject::PRESENT; LOG(info) << mID.getName() << " simulated in " << ((mROMode == o2::parameters::GRPObject::CONTINUOUS) ? "CONTINUOUS" : "TRIGGERED")