should use the total eff but not only seed

Author: ChunzhengLab

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

@@ -57,25 +57,25 @@ 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";
+  std::cout << "=====================\n";
   LOGP(info, "ALPIDE Vbb=0V response");
   mAlpSimResp0 = loadResponse(AlpideFile, "response0"); // Vbb=0V
   mAlpSimResp0->computeCentreFromData();
   mAlpSimResp0->print();
   LOGP(info, "Response Centre {}", mAlpSimResp0->getRespCentreDep());
-  std::cout<<"=====================\n";
+  std::cout << "=====================\n";
   LOGP(info, "ALPIDE Vbb=-3V response");
   mAlpSimResp1 = loadResponse(AlpideFile, "response1"); // Vbb=-3V
   mAlpSimResp1->computeCentreFromData();
   mAlpSimResp1->print();
   LOGP(info, "Response Centre {}", mAlpSimResp1->getRespCentreDep());
-  std::cout<<"=====================\n";
+  std::cout << "=====================\n";
   LOGP(info, "APTS response");
   mAptSimResp1 = loadResponse(AptsFile, "response1"); // APTS
   mAptSimResp1->computeCentreFromData();
   mAptSimResp1->print();
   LOGP(info, "Response Centre {}", mAptSimResp1->getRespCentreDep());
-  std::cout<<"=====================\n";
+  std::cout << "=====================\n";
 }
 
 std::vector<float> getCollectionSeediciencies(o2::its3::ChipSimResponse* resp,
@@ -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,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);
 
-  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;
+    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 = 0.f; return; }
+
+  if (totQ <= 0.f) {
+    mRespCentreDep = mDptMin;
+    return;
+  }
 
   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);
+  for (auto& b : bins) {
+    if (cumQ + b.dq < halfQ) {
+      cumQ += b.dq;
+      continue;
+    }
     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 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(A) > 1.e-12f)
-      x = (-B + std::sqrt(disc)) / (2.f * A);
-    else
+    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 = mDptMax;
 }