diff --git a/Detectors/Upgrades/ITS3/macros/test/CheckChipResponseFile.C b/Detectors/Upgrades/ITS3/macros/test/CheckChipResponseFile.C index 996a99d87ecbc..32d5bad87ce21 100644 --- a/Detectors/Upgrades/ITS3/macros/test/CheckChipResponseFile.C +++ b/Detectors/Upgrades/ITS3/macros/test/CheckChipResponseFile.C @@ -24,6 +24,7 @@ #define ENABLE_UPGRADES #include "ITSMFTSimulation/AlpideSimResponse.h" +#include "ITS3Simulation/ChipSimResponse.h" #include "ITS3Base/SegmentationMosaix.h" #include "fairlogger/Logger.h" @@ -34,21 +35,21 @@ 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; +std::unique_ptr mAlpSimResp0, mAlpSimResp1, mAptSimResp1; -o2::itsmft::AlpideSimResponse* loadResponse(const std::string& fileName, const std::string& respName) +std::unique_ptr 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) + auto base = f->Get(respName.c_str()); + if (!base) { std::cerr << respName << " not found in " << fileName << std::endl; - return resp; + return nullptr; + } + return std::make_unique(base); } void LoadRespFunc() @@ -56,39 +57,49 @@ 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"; + std::cout << "=====================\n"; + LOGP(info, "ALPIDE Vbb=0V response"); mAlpSimResp0 = loadResponse(AlpideFile, "response0"); // Vbb=0V - LOG(info) << "ALPIDE Vbb=0V response" << std::endl; + mAlpSimResp0->computeCentreFromData(); mAlpSimResp0->print(); + LOGP(info, "Response Centre {}", mAlpSimResp0->getRespCentreDep()); + std::cout << "=====================\n"; + LOGP(info, "ALPIDE Vbb=-3V response"); mAlpSimResp1 = loadResponse(AlpideFile, "response1"); // Vbb=-3V - LOG(info) << "ALPIDE Vbb=-3V response" << std::endl; + mAlpSimResp1->computeCentreFromData(); mAlpSimResp1->print(); + LOGP(info, "Response Centre {}", mAlpSimResp1->getRespCentreDep()); + std::cout << "=====================\n"; + LOGP(info, "APTS response"); mAptSimResp1 = loadResponse(AptsFile, "response1"); // APTS - LOG(info) << "APTS response" << std::endl; + mAptSimResp1->computeCentreFromData(); mAptSimResp1->print(); + LOGP(info, "Response Centre {}", mAptSimResp1->getRespCentreDep()); + std::cout << "=====================\n"; } -std::vector getCollectionSeediciencies(o2::itsmft::AlpideSimResponse* resp, +std::vector getCollectionSeediciencies(o2::its3::ChipSimResponse* 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, + um2cm(depth) + 1.e-9, flipRow, flipCol); seed.push_back(rspmat ? rspmat->getValue(2, 2) : 0.f); } return seed; } -std::vector getShareValues(o2::itsmft::AlpideSimResponse* resp, +std::vector getShareValues(o2::its3::ChipSimResponse* 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, + um2cm(depth) + 1.e-9, flipRow, flipCol); float s = 0; int npix = resp->getNPix(); @@ -103,14 +114,14 @@ std::vector getShareValues(o2::itsmft::AlpideSimResponse* resp, return share; } -std::vector getEffValues(o2::itsmft::AlpideSimResponse* resp, +std::vector getEffValues(o2::its3::ChipSimResponse* 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, + um2cm(depth) + 1.e-9, flipRow, flipCol); float s = 0; int npix = resp->getNPix(); @@ -129,13 +140,16 @@ void CheckChipResponseFile() LoadRespFunc(); LOG(info) << "Response function loaded" << std::endl; - std::vector vecDepth(50); - for (int i = 0; i < 50; ++i) - vecDepth[i] = i; + std::vector vecDepth; + int numPoints = 100; + for (int i = 0; i < numPoints; ++i) { + float value = -50 + i * (100.0f / (numPoints - 1)); + vecDepth.push_back(value); + } int colors[] = {kOrange + 7, kRed + 1, kAzure + 4}; struct RespInfo { - o2::itsmft::AlpideSimResponse* resp; + std::unique_ptr& resp; std::string title; int color; }; @@ -145,7 +159,7 @@ void CheckChipResponseFile() {mAlpSimResp1, "ALPIDE Vbb=-3V", colors[2]}}; TCanvas* c1 = new TCanvas("c1", "c1", 800, 600); - TH1* frame = c1->DrawFrame(-1, -0.049, 50, 1.049); + TH1* frame = c1->DrawFrame(-50, -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); @@ -154,11 +168,11 @@ void CheckChipResponseFile() 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); + auto seed = getCollectionSeediciencies(r.resp.get(), vecDepth); + auto shr = getShareValues(r.resp.get(), vecDepth); + auto all = getEffValues(r.resp.get(), vecDepth); - TGraph* grSeed = new TGraph(vecDepth.size(), vecDepth.data(), seed.data()); + auto 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); @@ -168,7 +182,7 @@ void CheckChipResponseFile() 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()); + auto grShare = new TGraph(vecDepth.size(), vecDepth.data(), shr.data()); grShare->SetLineColor(r.color); grShare->SetLineWidth(2); grShare->SetMarkerColor(r.color); @@ -177,7 +191,7 @@ void CheckChipResponseFile() 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()); + auto grEff = new TGraph(vecDepth.size(), vecDepth.data(), all.data()); grEff->SetLineColor(r.color); grEff->SetLineWidth(2); grEff->SetMarkerColor(r.color); diff --git a/Detectors/Upgrades/ITS3/simulation/src/ChipDigitsContainer.cxx b/Detectors/Upgrades/ITS3/simulation/src/ChipDigitsContainer.cxx index 102b15863683e..efe878536687d 100644 --- a/Detectors/Upgrades/ITS3/simulation/src/ChipDigitsContainer.cxx +++ b/Detectors/Upgrades/ITS3/simulation/src/ChipDigitsContainer.cxx @@ -31,12 +31,12 @@ void ChipDigitsContainer::addNoise(UInt_t rofMin, UInt_t rofMax, const o2::its3: int nel = 0; if (isIB()) { - // Inner barrel: use ITS3-specific noise interface with OB segmentation. - mean = params->getIBNoisePerPixel() * SegmentationOB::NPixels; + // Inner barrel: use ITS3-specific noise interface with IB segmentation. + mean = params->getIBNoisePerPixel() * SegmentationIB::NPixels; nel = static_cast(params->getIBChargeThreshold() * 1.1); } else { - // Outer barrel: use base class noise interface with IB segmentation. - mean = params->getNoisePerPixel() * SegmentationIB::NPixels; + // Outer barrel: use base class noise interface with OB segmentation. + mean = params->getNoisePerPixel() * SegmentationOB::NPixels; nel = static_cast(params->getChargeThreshold() * 1.1); } diff --git a/Detectors/Upgrades/ITS3/simulation/src/ChipSimResponse.cxx b/Detectors/Upgrades/ITS3/simulation/src/ChipSimResponse.cxx index 1c482983f0d0a..72b291fb0d653 100644 --- a/Detectors/Upgrades/ITS3/simulation/src/ChipSimResponse.cxx +++ b/Detectors/Upgrades/ITS3/simulation/src/ChipSimResponse.cxx @@ -25,38 +25,66 @@ void ChipSimResponse::initData(int tableNumber, std::string dataPath, const bool void ChipSimResponse::computeCentreFromData() { - std::vector zVec, qVec; const int npix = o2::itsmft::AlpideRespSimMat::getNPix(); + std::vector zVec, effVec; + zVec.reserve(mNBinDpt); + effVec.reserve(mNBinDpt); 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); + int rev = mNBinDpt - 1 - iz; + float z = mDptMin + iz / mStepInvDpt; + float sum = 0.f; + const auto& mat = mData[rev]; + for (int ix = 0; ix < npix; ++ix) { + for (int iy = 0; iy < npix; ++iy) { + sum += mat.getValue(ix, iy); + } + } + zVec.push_back(z); + effVec.push_back(sum); } - 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); - } + struct Bin { + float z0, z1, q0, q1, dq; + }; + std::vector bins; + bins.reserve(zVec.size() - 1); - float intQ = 0.f, intZQ = 0.f; + float totQ = 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; + float q0 = effVec[i], q1 = effVec[i + 1]; + float dq = 0.5f * (q0 + q1) * (z1 - z0); + bins.push_back({z0, z1, q0, q1, dq}); + totQ += dq; + } + + if (totQ <= 0.f) { + mRespCentreDep = mDptMin; + return; + } + + float halfQ = 0.5f * totQ; + float cumQ = 0.f; + for (auto& b : bins) { + if (cumQ + b.dq < halfQ) { + cumQ += b.dq; + continue; + } + float dz = b.z1 - b.z0; + float slope = (b.q1 - b.q0) / dz; + float disc = b.q0 * b.q0 - 2.f * slope * (cumQ - halfQ); + + float x; + if (disc >= 0.f && std::abs(slope) > 1e-6f) { + x = (-b.q0 + std::sqrt(disc)) / slope; + } else { + x = (halfQ - cumQ) / b.q0; + } + x = std::clamp(x, 0.f, dz); + mRespCentreDep = b.z0 + x; + return; } - mRespCentreDep = (intQ > 0.f) ? intZQ / intQ : 0.f; + mRespCentreDep = mDptMax; }