sPhenix code displayed by LXR master/​coresoftware/​offline/​packages/​CaloReco/​BEmcRecCEMC.cc	Source navigation Diff markup Identifier search General search
	Version: master Revert   Change

File indexing completed on 2025-12-17 09:19:56

 #include "BEmcRecCEMC.h"
 
 #include "BEmcCluster.h"
 #include "BEmcProfile.h"
 
 #include <cmath>
 #include <iostream>
 #include <Math/Vector3D.h>   // for TVector3 used in calculateIncidence
 #include <algorithm>    // for std::max
 using ROOT::Math::XYZVector;
 
 BEmcRecCEMC::BEmcRecCEMC()
 {
   Name("BEmcRecCEMC");
   SetCylindricalGeometry();
 }
 
 void BEmcRecCEMC::LoadProfile(const std::string& fname)
 {
   //  std::cout << "Infor from BEmcRecCEMC::LoadProfile(): no external file used for shower profile evaluation in CEMC" << std::endl;
   _emcprof = new BEmcProfile(fname);
 }
 
 void BEmcRecCEMC::GetImpactThetaPhi(float xg, float yg, float zg, float& theta, float& phi)
 {
   theta = 0;
   phi = 0;
 
   //  float theta = atan(sqrt(xg*xg + yg*yg)/fabs(zg-fVz));
   float rg = std::sqrt((xg * xg) + (yg * yg));
   float theta_twr;
   if (std::fabs(zg) <= 15)
   {
     theta_twr = 0;
   }
   else if (zg > 15)
   {
     theta_twr = std::atan2(zg - 15, rg);
   }
   else
   {
     theta_twr = std::atan2(zg + 15, rg);
   }
   float theta_tr = std::atan2(zg - fVz, rg);
   theta = std::fabs(theta_tr - theta_twr);
   //  phi = atan2(yg,xg);
 }
 
 /*
 float BEmcRecCEMC::GetProb(vector<EmcModule> HitList, float ecl, float xg, float yg, float zg, float& chi2, int& ndf)
 {
   chi2 = 0;
   ndf = 0;
   float prob = -1;
 
   //  float theta = atan(sqrt(xg*xg + yg*yg)/fabs(zg-fVz));
   float rg = sqrt(xg * xg + yg * yg);
   float theta_twr;
   if (fabs(zg) <= 15)
     theta_twr = 0;
   else if (zg > 15)
     theta_twr = atan2(zg - 15, rg);
   else
     theta_twr = atan2(zg + 15, rg);
   float theta_tr = atan2(zg - fVz, rg);
   float theta = fabs(theta_tr - theta_twr);
 
   float phi = atan2(yg, xg);
   if (_emcprof != nullptr) prob = _emcprof->GetProb(&HitList, fNx, ecl, theta, phi);
 
   return prob;
 }
 */
 /*
 float BEmcRecCEMC::GetProb(vector<EmcModule> HitList, float et, float xg, float yg, float zg, float& chi2, int& ndf)
 // et, xg, yg, zg not used here
 {
   const float thresh = 0.01;
   const int DXY = 3;  // 2 is for 5x5 matrix; 3 for 7x7 matrix
   const int Nmax = 1000;
   float ee[Nmax];
   int iyy[Nmax];
   int izz[Nmax];
 
   int ich;
   vector<EmcModule>::iterator ph = HitList.begin();
 
   chi2 = 0;
   ndf = 0;
 
   int nn = 0;
 
   while (ph != HitList.end())
   {
     ee[nn] = ph->amp;
     if (ee[nn] > thresh)
     {
       ich = ph->ich;
       izz[nn] = ich % fNx;
       iyy[nn] = ich / fNx;
       nn++;
       if (nn >= Nmax)
       {
       std::cout << "BEmcRec::GetProb: Cluster size is too big. Skipping the rest of the towers" << std::endl;
         break;
       }
     }  // if( ee[nn]
     ++ph;
   }  // while( ph
 
   if (nn <= 0) return -1;
 
   int iy0 = -1, iz0 = -1;
   float emax = 0;
 
   for (int i = 0; i < nn; i++)
   {
     if (ee[i] > emax)
     {
       emax = ee[i];
       iy0 = iyy[i];
       iz0 = izz[i];
     }
   }
 
   if (emax <= 0) return -1;
 
   int id;
   float etot = 0;
   float sz = 0;
   float sy = 0;
 
   for (int idz = -DXY; idz <= DXY; idz++)
   {
     for (int idy = -DXY; idy <= DXY; idy++)
     {
       id = GetTowerID(iy0 + idy, iz0 + idz, nn, iyy, izz, ee);
       if (id >= 0)
       {
         etot += ee[id];
         sz += ee[id] * (iz0 + idz);
         sy += ee[id] * (iy0 + idy);
       }
     }
   }
   float zcg = sz / etot;  // Here cg allowed to be out of range
   float ycg = sy / etot;
   int iz0cg = int(zcg + 0.5);
   int iy0cg = int(ycg + 0.5);
   float ddz = fabs(zcg - iz0cg);
   float ddy = fabs(ycg - iy0cg);
 
   int isz = 1;
   if (zcg - iz0cg < 0) isz = -1;
   int isy = 1;
   if (ycg - iy0cg < 0) isy = -1;
 
   // 4 central towers: 43
   //                   12
   // Tower 1 - central one
   float e1, e2, e3, e4;
   e1 = e2 = e3 = e4 = 0;
   id = GetTowerID(iy0cg, iz0cg, nn, iyy, izz, ee);
   if (id >= 0) e1 = ee[id];
   id = GetTowerID(iy0cg, iz0cg + isz, nn, iyy, izz, ee);
   if (id >= 0) e2 = ee[id];
   id = GetTowerID(iy0cg + isy, iz0cg + isz, nn, iyy, izz, ee);
   if (id >= 0) e3 = ee[id];
   id = GetTowerID(iy0cg + isy, iz0cg, nn, iyy, izz, ee);
   if (id >= 0) e4 = ee[id];
 
   float e1t = (e1 + e2 + e3 + e4) / etot;
   float e2t = (e1 + e2 - e3 - e4) / etot;
   float e3t = (e1 - e2 - e3 + e4) / etot;
   float e4t = (e3) / etot;
   //  float e5t = (e2+e4)/etot;
 
   float rr = sqrt((0.5 - ddz) * (0.5 - ddz) + (0.5 - ddy) * (0.5 - ddy));
 
   float c1, c2, c11;
 
   float logE = log(etot);
 
   // e1 energy is the most effective for PID if properly tuned !
   // Discrimination power is very sensitive to paramter c1: the bigger it is
   // the better discrimination;
   c1 = 0.95;
   c2 = 0.0066364 * logE + 0.00466667;
   if (c2 < 0) c2 = 0;
   float e1p = c1 - c2 * rr * rr;
   c1 = 0.034 - 0.01523 * logE + 0.0029 * logE * logE;
   float err1 = c1;
 
   // For e2
   c1 = 0.00844086 + 0.00645359 * logE - 0.00119381 * logE * logE;
   if (etot > 15) c1 = 0.00844086 + 0.00645359 * log(15.) - 0.00119381 * log(15.) * log(15.);  // Const at etot>15GeV
   if (c1 < 0) c1 = 0;
   c2 = 3.9;                                                      // Fixed
   float e2p = sqrt(c1 + 0.25 * c2) - sqrt(c1 + c2 * ddy * ddy);  // =0 at ddy=0.5
 
   c1 = 0.0212333 + 0.0420473 / etot;
   c2 = 0.090;  // Fixed
   float err2 = c1 + c2 * ddy;
   if (ddy > 0.3) err2 = c1 + c2 * 0.3;  // Const at ddy>0.3
 
   // For e3
   c1 = 0.0107857 + 0.0056801 * logE - 0.000892016 * logE * logE;
   if (etot > 15) c1 = 0.0107857 + 0.0056801 * log(15.) - 0.000892016 * log(15.) * log(15.);  // Const at etot>15GeV
   if (c1 < 0) c1 = 0;
   c2 = 3.9;                                                      // Fixed
   float e3p = sqrt(c1 + 0.25 * c2) - sqrt(c1 + c2 * ddz * ddz);  // =0 at ddz=0.5
 
   //  c1 = 0.0200 + 0.042/etot;
   c1 = 0.0167 + 0.058 / etot;
   c2 = 0.090;  // Fixed
   float err3 = c1 + c2 * ddz;
   if (ddz > 0.3) err3 = c1 + c2 * 0.3;  // Const at ddz>0.3
 
   // For e4
   float e4p = 0.25 - 0.668 * rr + 0.460 * rr * rr;
   c11 = 0.171958 + 0.0142421 * logE - 0.00214827 * logE * logE;
   //  c11 = 0.171085 + 0.0156215*logE - -0.0025809*logE*logE;
   float err4 = 0.102 - 1.43 * c11 * rr + c11 * rr * rr;  // Min is set to x=1.43/2.
   err4 *= 1.1;
 
   chi2 = 0.;
   chi2 += (e1p - e1t) * (e1p - e1t) / err1 / err1;
   chi2 += (e2p - e2t) * (e2p - e2t) / err2 / err2;
   chi2 += (e3p - e3t) * (e3p - e3t) / err3 / err3;
   chi2 += (e4p - e4t) * (e4p - e4t) / err4 / err4;
   ndf = 4;
 
   //  chi2 /= 1.1;
   float prob = TMath::Prob(chi2, ndf);
 
   return prob;
 }
 */
 
 bool BEmcRecCEMC::calculateIncidence(float x, float y,
                                      float& a_phi_sgn, float& a_eta_sgn)
 {
   // ────────────────────────────────────────────────────────────────────────────
   // Step 1) Owner tower in "tower units"
   // ────────────────────────────────────────────────────────────────────────────
   const int ix = EmcCluster::lowint(x + 0.5F);
   const int iy = EmcCluster::lowint(y + 0.5F);
 
   TowerGeom g{};
   if (!GetTowerGeometry(ix, iy, g))
   {
     return false;  // need detailed geometry
   }
 
   // Sub-cell offsets δ ∈ (−0.5, +0.5]
   const float dx = x - ix;
   const float dy = y - iy;
 
   // ────────────────────────────────────────────────────────────────────────────
   // Step 2) Front-surface point F in cm:
   //         F = C + δx * eφ_geom + δy * eη_geom
   // ────────────────────────────────────────────────────────────────────────────
   const XYZVector C(g.Xcenter, g.Ycenter, g.Zcenter);
   const XYZVector ephi_geom(g.dX[0], g.dY[0], g.dZ[0]);
   const XYZVector eeta_geom(g.dX[1], g.dY[1], g.dZ[1]);
 
   if (std::sqrt(ephi_geom.Mag2()) < 1e-9 || std::sqrt(eeta_geom.Mag2()) < 1e-9)
   {
     return false;
   }
 
   const XYZVector F = C + dx * ephi_geom + dy * eeta_geom;
 
   // ────────────────────────────────────────────────────────────────────────────
   // Step 3) Unit ray from vertex to the front surface
   // ────────────────────────────────────────────────────────────────────────────
   const XYZVector V(0., 0., fVz);
   XYZVector pF = F - V;
   const double pFmag = std::sqrt(pF.Mag2());
   if (!(pFmag > 0.0))
   {
     return false;
   }
   pF *= (1.0 / pFmag);
 
   // Helper: from a basis {un,uphi,ueta} and ray pF, compute projections
   // and signed incidence angles (no cosines returned).
   auto basis_to_incidence =
     [&](const XYZVector& un_b, const XYZVector& uphi_b, const XYZVector& ueta_b,
         double& pn_b, double& pph_b, double& pet_b,
         double& aphi_b, double& aeta_b) -> bool
   {
     pn_b  = pF.Dot(un_b);
     pph_b = pF.Dot(uphi_b);
     pet_b = pF.Dot(ueta_b);
 
     const double apn_b      = std::max(1e-12, std::fabs(pn_b));
     const double aphi_mag_b = std::atan2(std::fabs(pph_b), apn_b);
     const double aeta_mag_b = std::atan2(std::fabs(pet_b), apn_b);
 
     const double sgn_phi_b = (pph_b >= 0.0) ? +1.0 : -1.0;
     const double sgn_eta_b = (pet_b >= 0.0) ? +1.0 : -1.0;
 
     aphi_b = sgn_phi_b * aphi_mag_b;
     aeta_b = sgn_eta_b * aeta_mag_b;
 
     return (std::isfinite(aphi_b) && std::isfinite(aeta_b));
   };
 
   // ────────────────────────────────────────────────────────────────────────────
   // Step 4) MECHANICAL frame from RawTowerGeomv5 rotations
   // ────────────────────────────────────────────────────────────────────────────
   XYZVector un_mech;
   XYZVector uphi_mech;
   XYZVector ueta_mech;
   {
     const double rx = g.rotX;
     const double ry = g.rotY;
     const double rz = g.rotZ;
 
     const double cx = std::cos(rx);
     const double sx = std::sin(rx);
     const double cy = std::cos(ry);
     const double sy = std::sin(ry);
     const double cz = std::cos(rz);
     const double sz = std::sin(rz);
 
     auto applyRot = [&](const XYZVector& v) -> XYZVector
     {
       // Apply Rz * Ry * Rx to local vector v (column-vector convention).
       double vx = v.X();
       double vy = v.Y();
       double vz = v.Z();
 
       // Rx
       const double vx1 = vx;
       const double vy1 =  cx * vy - sx * vz;
       const double vz1 =  sx * vy + cx * vz;
 
       // Ry
       const double vx2 =  cy * vx1 + sy * vz1;
       const double vy2 =  vy1;
       const double vz2 = -sy * vx1 + cy * vz1;
 
       // Rz
       const double vx3 =  cz * vx2 - sz * vy2;
       const double vy3 =  sz * vx2 + cz * vy2;
       const double vz3 =  vz2;
 
       return XYZVector(vx3, vy3, vz3);
     };
 
     // Local axes: ẑ = bar axis, x̂/ŷ = transverse directions
     XYZVector u_axis = applyRot(XYZVector(0., 0., 1.)); // tower "z" axis
     const XYZVector u_x    = applyRot(XYZVector(1., 0., 0.)); // tower "x" axis
     const XYZVector u_y    = applyRot(XYZVector(0., 1., 0.)); // tower "y" axis
 
     // Normalize axis and enforce outward direction
     const double amag = std::sqrt(u_axis.Mag2());
     if (!(amag > 0.0))
     {
       return false;
     }
 
     u_axis *= (1.0 / amag);
     if (u_axis.Dot(C) < 0.0)
     {
       u_axis = -u_axis;
     }
 
     un_mech = u_axis;
 
     // φ̂: project rotated local x̂ into plane ⟂ axis; normalize
     uphi_mech = u_x - un_mech * u_x.Dot(un_mech);
     const double uphim = std::sqrt(uphi_mech.Mag2());
     if (!(uphim > 0.0))
     {
       return false;
     }
     uphi_mech *= (1.0 / uphim);
 
     // η̂: start from rotated local ŷ, orthogonalize to {un,uphi}; normalize
     ueta_mech = u_y
               - un_mech   * u_y.Dot(un_mech)
               - uphi_mech * u_y.Dot(uphi_mech);
     const double uetam = std::sqrt(ueta_mech.Mag2());
     if (!(uetam > 0.0))
     {
       return false;
     }
     ueta_mech *= (1.0 / uetam);
 
     // Ensure right-handed triad
     if (un_mech.Cross(uphi_mech).Dot(ueta_mech) < 0.0)
     {
       ueta_mech = -ueta_mech;
     }
 
     // Align mechanical φ-axis sign with geometric φ tangent if available
     XYZVector uphi_ref = ephi_geom;
     const double refMag = std::sqrt(uphi_ref.Mag2());
     if (refMag > 0.0)
     {
       uphi_ref *= (1.0 / refMag);
       if (uphi_mech.Dot(uphi_ref) < 0.0)
       {
         uphi_mech = -uphi_mech;
         ueta_mech = -ueta_mech;
       }
     }
   }
 
   // ────────────────────────────────────────────────────────────────────────────
   // Step 5) Compute incidence in the mechanical frame
   // ────────────────────────────────────────────────────────────────────────────
   double pn   = 0.0;
   double pph  = 0.0;
   double pet  = 0.0;
   double aphi = 0.0;
   double aeta = 0.0;
 
   if (!basis_to_incidence(un_mech, uphi_mech, ueta_mech,
                           pn, pph, pet,
                           aphi, aeta))
   {
     return false;
   }
 
   a_phi_sgn = static_cast<float>(aphi);
   a_eta_sgn = static_cast<float>(aeta);
 
   m_lastAlphaPhiSigned = a_phi_sgn;
   m_lastAlphaEtaSigned = a_eta_sgn;
 
   const bool ok =
       (std::isfinite(a_phi_sgn) && std::isfinite(a_eta_sgn));
 
   return ok;
 }
 
 void BEmcRecCEMC::CorrectShowerDepth(int ix, int iy, float E, float xA, float yA, float zA, float& xC, float& yC, float& zC)
 {
   if (!m_UseCorrectShowerDepth)
   {
     xC = xA;
     yC = yA;
     zC = zA;
     return;
   }
 
   float logE = log(0.1);
   if (E > 0.1)
   {
     logE = std::log(E);
   }
 
   float phi = 0;
   // Rotate by phi (towers are tilted by a fixed angle in phi by ~9 deg?)
   // Just tuned from sim data
   if (m_UseDetailedGeometry)
   {
     phi = 0.0033 - 0.0010 * logE;
   }
   else
   {
     phi = 0.0016 - 0.0010 * logE;
   }
   xC = xA * std::cos(phi) - yA * std::sin(phi);
   yC = xA * std::sin(phi) + yA * std::cos(phi);
 
   // Correction in z
   float rA = std::sqrt((xA * xA) + (yA * yA));
   float theta_twr;
   if (m_UseDetailedGeometry)
   {
     // Read the angle right from the detailed RawTowerGeom objects
     TowerGeom geom0;
     GetTowerGeometry(ix, iy, geom0);
     theta_twr = M_PI / 2 + geom0.rotX;
   }
   else
   {
     // Use the approximate default tower angles array.
     theta_twr = M_PI / 2 - angles[iy];
   }
 
   float theta_tr = std::atan2(zA - fVz, rA);
 
   // Shower CG in long. direction
   // Different tuning for the approximate and detailed geometry
   float L = 0;
   if (m_UseDetailedGeometry)
   {
     L = -2.67787 + 0.924138 * logE;
   }
   else
   {
     L = -1.79968 + 0.837322 * logE;
   }
   float dz = L * std::sin(theta_tr - theta_twr) / std::cos(theta_twr);
 
   if (!m_UseDetailedGeometry)
   {
     // Not strictly speaking a "shower depth correction" but rather a baseline correction
     // The factor 0.10 accounts for the fact that the approximate geometry
     // is projected at 93.5cm, which is roughly 90 % of the actual average EMCal radius (~ 105 cm)
     dz -= fVz * 0.10;
   }
 
   zC = zA - dz;
 
   // zvtx-dependent pseudorapidity-correction
   // At large zvtx (> 20 cm), the sawtooth shape of the EMCal is no longer projective wrt collision point,
   // it introduces a bias which can be approximately corrected with the linear formula below
   if (m_UseDetailedGeometry && std::abs(fVz) > 20)
   {
     float eta = std::asinh((zC - fVz) / rA);
     if (fVz > 0)
     {
       eta -= slopes_[iy] * (fVz - 20);
     }
     else
     {
       eta -= slopes_[fNy - iy - 1] * (fVz + 20);
     }
     zC = fVz + rA * std::sinh(eta);
   }
 
   return;
 }
 
 void BEmcRecCEMC::CorrectEnergy(float Energy, float /*x*/, float /*y*/,
                                 float& Ecorr)
 {
   // Corrects the EM Shower Energy for attenuation in fibers and
   // long energy leakage
   //
   // (x,y) - impact position (cm) in Sector frame
   /*
   float sinT;
   float att, leak, corr;
   const float leakPar = 0.0033; // parameter from fit
   const float attPar = 120; // Attenuation in module (cm)
   const float X0 = 2; // radiation length (cm)
 
   *Ecorr = Energy;
   if( Energy < 0.01 ) return;
 
   GetImpactAngle(x, y, &sinT); // sinT not used so far
   leak = 2-sqrt(1+leakPar*log(1+Energy)*log(1+Energy));
   att = exp(log(Energy)*X0/attPar);
   corr = leak*att;
   *Ecorr = Energy/corr;
   */
   Ecorr = Energy;
 }
 
 void BEmcRecCEMC::CorrectECore(float Ecore, float /*x*/, float /*y*/, float& Ecorr)
 {
   // Corrects the EM Shower Core Energy for attenuation in fibers,
   // long energy leakage and angle dependance
   //
   // (x,y) - impact position (cm) in Sector frame
 
   const float c0 = 0.950;  // For no threshold
   Ecorr = Ecore / c0;
 }
 
 void BEmcRecCEMC::CorrectPosition(float Energy, float x, float y,
                                   float& xc, float& yc)
 {
   // Corrects the Shower Center of Gravity for the systematic shift due to
   // limited tower size
   //
   // Everything here is in tower units.
   // (x,y) - CG position, (xc,yc) - corrected position
 
   float bx;
   float by;
   float x0;
   float y0;
   int ix0;
   int iy0;
 
   if (!m_UseCorrectPosition)
   {
     xc = x;
     yc = y;
     return;
   }
 
   if (Energy < 0.01)
   {
     return;
   }
 
   bx = 0.15;
   by = 0.15;
 
   x0 = x;
   ix0 = EmcCluster::lowint(x0 + 0.5);
 
   if (std::abs(x0 - ix0) <= 0.5)
   {
     x0 = ix0 + bx * asinh(2. * (x0 - ix0) * sinh(0.5 / bx));
   }
   else
   {
     x0 = x;
     std::cout << "????? Something wrong in BEmcRecCEMC::CorrectPosition: x = "
               << x << " dx = " << x0 - ix0 << std::endl;
   }
 
   // Correct for phi bias within module of 8 towers
 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
   int ix8 = int(x + 0.5) / 8; // that is hokey - suggest lroundf(x)
   float x8 = x + 0.5 - (ix8 * 8) - 4;  // from -4 to +4
   float dx = 0;
   if (m_UseDetailedGeometry)
   {
     // Don't know why there is a different factor for each tower of the sector
     // Just tuned from MC
 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
     int local_ix8 = int(x+0.5) - ix8 * 8; // that is hokey - suggest lroundf(x)
     dx = factor_[local_ix8] * x8 / 4.;
   }
   else
   {
     dx = 0.10 * x8 / 4.;
     if (std::fabs(x8) > 3.3)
     {
       dx = 0;  // Don't correct near the module edge
     }
   }
 
   xc = x0 - dx;
   while (xc < -0.5)
   {
     xc += float(fNx);
   }
   while (xc >= fNx - 0.5)
   {
     xc -= float(fNx);
   }
 
   y0 = y;
   iy0 = EmcCluster::lowint(y0 + 0.5);
 
   if (std::abs(y0 - iy0) <= 0.5)
   {
     y0 = iy0 + by * std::asinh(2. * (y0 - iy0) * sinh(0.5 / by));
   }
   else
   {
     y0 = y;
     std::cout << "????? Something wrong in BEmcRecCEMC::CorrectPosition: y = "
               << y << "dy = " << y0 - iy0 << std::endl;
   }
   yc = y0;
 }

This page was automatically generated by the 2.3.7 LXR engine. The LXR team