should use the total eff but not only seed

Author: ChunzhengLab

Detectors/Upgrades/ITS3/macros/test/CheckChipResponseFile.C

@@ -85,7 +85,7 @@ std::vector<float> getCollectionSeediciencies(o2::its3::ChipSimResponse* resp,
   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);
   }
@@ -99,7 +99,7 @@ std::vector<float> getShareValues(o2::its3::ChipSimResponse* resp,
   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();
@@ -121,7 +121,7 @@ std::vector<float> getEffValues(o2::its3::ChipSimResponse* resp,
   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();
@@ -140,9 +140,12 @@ 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 {
@@ -156,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);

Detectors/Upgrades/ITS3/simulation/src/ChipSimResponse.cxx

@@ -25,61 +25,62 @@ 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();
+    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);
-  }
+    for (int iz = 0; iz < mNBinDpt; ++iz) {
+        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);
 
-  struct BinInfo { float z0, z1, q0, q1, dq; };
-  std::vector<BinInfo> bins;
-  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;
-    float dq = 0.5f * (q0 + q1) * dz;
-    bins.push_back({z0, z1, q0, q1, dq});
-    totQ += dq;
-  }
-  if (totQ <= 0.f) { mRespCentreDep = 0.f; return; }
+    float totQ = 0.f;
+    for (size_t i = 0; i + 1 < zVec.size(); ++i) {
+        float z0 = zVec[i], z1 = zVec[i + 1];
+        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;
+    }
 
-  float halfQ = 0.5f * totQ;
-  float cumQ = 0.f;
-  for (const auto& b : bins) {
-    if (cumQ + b.dq < halfQ) { cumQ += b.dq; continue; }
-    float qSlope = (b.q1 - b.q0) / (b.z1 - b.z0);
-    float dz = b.z1 - b.z0;
-    float A = qSlope * 0.5f;
-    float B = b.q0;
-    float C = cumQ - halfQ;
-    float disc = B * B - 4.f * A * C;
-    float x;
-    if (disc >= 0.f && std::abs(A) > 1.e-12f)
-      x = (-B + std::sqrt(disc)) / (2.f * A);
-    else
-      x = (halfQ - cumQ) / b.q0;
-    x = std::clamp(x, 0.f, dz);
-    mRespCentreDep = b.z0 + x;
-    return;
-  }
-  mRespCentreDep = mDptMax;
+    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-12f) {
+            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 = mDptMax;
 }