@@ -25,62 +25,64 @@ void ChipSimResponse::initData(int tableNumber, std::string dataPath, const bool
2525
2626void ChipSimResponse::computeCentreFromData ()
2727{
28- const int npix = o2::itsmft::AlpideRespSimMat::getNPix ();
29- std::vector<float > zVec, effVec;
30- zVec.reserve (mNBinDpt );
31- effVec.reserve (mNBinDpt );
28+ const int npix = o2::itsmft::AlpideRespSimMat::getNPix ();
29+ std::vector<float > zVec, effVec;
30+ zVec.reserve (mNBinDpt );
31+ effVec.reserve (mNBinDpt );
3232
33- for (int iz = 0 ; iz < mNBinDpt ; ++iz) {
34- int rev = mNBinDpt - 1 - iz;
35- float z = mDptMin + iz / mStepInvDpt ;
36- float sum = 0 .f ;
37- const auto & mat = mData [rev];
38- for (int ix = 0 ; ix < npix; ++ix)
39- for (int iy = 0 ; iy < npix; ++iy)
40- sum += mat.getValue (ix, iy);
41- zVec.push_back (z);
42- effVec.push_back (sum);
43- }
33+ for (int iz = 0 ; iz < mNBinDpt ; ++iz) {
34+ int rev = mNBinDpt - 1 - iz;
35+ float z = mDptMin + iz / mStepInvDpt ;
36+ float sum = 0 .f ;
37+ const auto & mat = mData [rev];
38+ for (int ix = 0 ; ix < npix; ++ix)
39+ for (int iy = 0 ; iy < npix; ++iy)
40+ sum += mat.getValue (ix, iy);
41+ zVec.push_back (z);
42+ effVec.push_back (sum);
43+ }
4444
45- struct Bin { float z0, z1, q0, q1, dq; };
46- std::vector<Bin> bins;
47- bins.reserve (zVec.size () - 1 );
45+ struct Bin {
46+ float z0, z1, q0, q1, dq;
47+ };
48+ std::vector<Bin> bins;
49+ bins.reserve (zVec.size () - 1 );
4850
49- float totQ = 0 .f ;
50- for (size_t i = 0 ; i + 1 < zVec.size (); ++i) {
51- float z0 = zVec[i], z1 = zVec[i + 1 ];
52- float q0 = effVec[i], q1 = effVec[i + 1 ];
53- float dq = 0 .5f * (q0 + q1) * (z1 - z0);
54- bins.push_back ({z0, z1, q0, q1, dq});
55- totQ += dq;
56- }
51+ float totQ = 0 .f ;
52+ for (size_t i = 0 ; i + 1 < zVec.size (); ++i) {
53+ float z0 = zVec[i], z1 = zVec[i + 1 ];
54+ float q0 = effVec[i], q1 = effVec[i + 1 ];
55+ float dq = 0 .5f * (q0 + q1) * (z1 - z0);
56+ bins.push_back ({z0, z1, q0, q1, dq});
57+ totQ += dq;
58+ }
5759
58- if (totQ <= 0 .f ) {
59- mRespCentreDep = mDptMin ;
60- return ;
61- }
60+ if (totQ <= 0 .f ) {
61+ mRespCentreDep = mDptMin ;
62+ return ;
63+ }
6264
63- float halfQ = 0 .5f * totQ;
64- float cumQ = 0 .f ;
65- for (auto & b : bins) {
66- if (cumQ + b.dq < halfQ) {
67- cumQ += b.dq ;
68- continue ;
69- }
70- float dz = b.z1 - b.z0 ;
71- float slope = (b.q1 - b.q0 ) / dz;
72- float disc = b.q0 * b.q0 - 2 .f * slope * (cumQ - halfQ);
65+ float halfQ = 0 .5f * totQ;
66+ float cumQ = 0 .f ;
67+ for (auto & b : bins) {
68+ if (cumQ + b.dq < halfQ) {
69+ cumQ += b.dq ;
70+ continue ;
71+ }
72+ float dz = b.z1 - b.z0 ;
73+ float slope = (b.q1 - b.q0 ) / dz;
74+ float disc = b.q0 * b.q0 - 2 .f * slope * (cumQ - halfQ);
7375
74- float x;
75- if (disc >= 0 .f && std::abs (slope) > 1e-12f ) {
76- x = (-b.q0 + std::sqrt (disc)) / slope;
77- } else {
78- x = (halfQ - cumQ) / b.q0 ;
79- }
80- x = std::clamp (x, 0 .f , dz);
81- mRespCentreDep = b.z0 + x;
82- return ;
76+ float x;
77+ if (disc >= 0 .f && std::abs (slope) > 1e-12f ) {
78+ x = (-b.q0 + std::sqrt (disc)) / slope;
79+ } else {
80+ x = (halfQ - cumQ) / b.q0 ;
8381 }
82+ x = std::clamp (x, 0 .f , dz);
83+ mRespCentreDep = b.z0 + x;
84+ return ;
85+ }
8486
85- mRespCentreDep = mDptMax ;
87+ mRespCentreDep = mDptMax ;
8688}
0 commit comments