ITS3: move the energy deposition wrt. centre of response (#14415)

ChunzhengLab · web-flow · commit 41279c2f65c0 · 2025-06-24T18:29:12.000+02:00
* move the energy deposition wrt. centre of response

* fix the order of IB/OB for noise seeding

* print out the resp centre in the CheckChipResponseFile

* find the real centre

* should use the total eff but not only seed
diff --git 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,61 +35,71 @@ 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<o2::its3::ChipSimResponse> mAlpSimResp0, mAlpSimResp1, mAptSimResp1;
 
-o2::itsmft::AlpideSimResponse* loadResponse(const std::string& fileName, const std::string& respName)
+std::unique_ptr<o2::its3::ChipSimResponse> 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<o2::itsmft::AlpideSimResponse>(respName.c_str());
+  if (!base) {
     std::cerr << respName << " not found in " << fileName << std::endl;
-  return resp;
+    return nullptr;
+  }
+  return std::make_unique<o2::its3::ChipSimResponse>(base);
 }
 
 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<float> getCollectionSeediciencies(o2::itsmft::AlpideSimResponse* resp,
+std::vector<float> getCollectionSeediciencies(o2::its3::ChipSimResponse* resp,
                                               const std::vector<float>& depths)
 {
   std::vector<float> 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<float> getShareValues(o2::itsmft::AlpideSimResponse* resp,
+std::vector<float> getShareValues(o2::its3::ChipSimResponse* resp,
                                   const std::vector<float>& depths)
 {
   std::vector<float> 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<float> getShareValues(o2::itsmft::AlpideSimResponse* resp,
   return share;
 }
 
-std::vector<float> getEffValues(o2::itsmft::AlpideSimResponse* resp,
+std::vector<float> getEffValues(o2::its3::ChipSimResponse* resp,
                                 const std::vector<float>& depths)
 {
   std::vector<float> 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<float> vecDepth(50);
-  for (int i = 0; i < 50; ++i)
-    vecDepth[i] = i;
+  std::vector<float> 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<o2::its3::ChipSimResponse>& 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
@@ -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<int>(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<int>(params->getChargeThreshold() * 1.1);
   }
 
diff --git 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<float> zVec, qVec;
   const int npix = o2::itsmft::AlpideRespSimMat::getNPix();
+  std::vector<float> 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<std::pair<float, float>> 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<Bin> 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;
 }
