packages/CaloReco/BEmcRecCEMC.cc

0001 #include "BEmcRecCEMC.h"
0002
0003 #include "BEmcCluster.h"
0004 #include "BEmcProfile.h"
0005
0006 #include <cmath>
0007 #include <iostream>
0008
0009 BEmcRecCEMC::BEmcRecCEMC()
0010 //  : _emcprof(nullptr)
0011 {
0012   Name("BEmcRecCEMC");
0013   SetCylindricalGeometry();
0014 }
0015
0016 void BEmcRecCEMC::LoadProfile(const std::string& fname)
0017 {
0018   //  std::cout << "Infor from BEmcRecCEMC::LoadProfile(): no external file used for shower profile evaluation in CEMC" << std::endl;
0019   _emcprof = new BEmcProfile(fname);
0020 }
0021
0022 void BEmcRecCEMC::GetImpactThetaPhi(float xg, float yg, float zg, float& theta, float& phi)
0023 {
0024   theta = 0;
0025   phi = 0;
0026
0027   //  float theta = atan(sqrt(xg*xg + yg*yg)/fabs(zg-fVz));
0028   float rg = std::sqrt((xg * xg) + (yg * yg));
0029   float theta_twr;
0030   if (std::fabs(zg) <= 15)
0031   {
0032     theta_twr = 0;
0033   }
0034   else if (zg > 15)
0035   {
0036     theta_twr = std::atan2(zg - 15, rg);
0037   }
0038   else
0039   {
0040     theta_twr = std::atan2(zg + 15, rg);
0041   }
0042   float theta_tr = std::atan2(zg - fVz, rg);
0043   theta = std::fabs(theta_tr - theta_twr);
0044   //  phi = atan2(yg,xg);
0045 }
0046
0047 /*
0048 float BEmcRecCEMC::GetProb(vector<EmcModule> HitList, float ecl, float xg, float yg, float zg, float& chi2, int& ndf)
0049 {
0050   chi2 = 0;
0051   ndf = 0;
0052   float prob = -1;
0053
0054   //  float theta = atan(sqrt(xg*xg + yg*yg)/fabs(zg-fVz));
0055   float rg = sqrt(xg * xg + yg * yg);
0056   float theta_twr;
0057   if (fabs(zg) <= 15)
0058     theta_twr = 0;
0059   else if (zg > 15)
0060     theta_twr = atan2(zg - 15, rg);
0061   else
0062     theta_twr = atan2(zg + 15, rg);
0063   float theta_tr = atan2(zg - fVz, rg);
0064   float theta = fabs(theta_tr - theta_twr);
0065
0066   float phi = atan2(yg, xg);
0067   if (_emcprof != nullptr) prob = _emcprof->GetProb(&HitList, fNx, ecl, theta, phi);
0068
0069   return prob;
0070 }
0071 */
0072 /*
0073 float BEmcRecCEMC::GetProb(vector<EmcModule> HitList, float et, float xg, float yg, float zg, float& chi2, int& ndf)
0074 // et, xg, yg, zg not used here
0075 {
0076   const float thresh = 0.01;
0077   const int DXY = 3;  // 2 is for 5x5 matrix; 3 for 7x7 matrix
0078   const int Nmax = 1000;
0079   float ee[Nmax];
0080   int iyy[Nmax];
0081   int izz[Nmax];
0082
0083   int ich;
0084   vector<EmcModule>::iterator ph = HitList.begin();
0085
0086   chi2 = 0;
0087   ndf = 0;
0088
0089   int nn = 0;
0090
0091   while (ph != HitList.end())
0092   {
0093     ee[nn] = ph->amp;
0094     if (ee[nn] > thresh)
0095     {
0096       ich = ph->ich;
0097       izz[nn] = ich % fNx;
0098       iyy[nn] = ich / fNx;
0099       nn++;
0100       if (nn >= Nmax)
0101       {
0102       std::cout << "BEmcRec::GetProb: Cluster size is too big. Skipping the rest of the towers" << std::endl;
0103         break;
0104       }
0105     }  // if( ee[nn]
0106     ++ph;
0107   }  // while( ph
0108
0109   if (nn <= 0) return -1;
0110
0111   int iy0 = -1, iz0 = -1;
0112   float emax = 0;
0113
0114   for (int i = 0; i < nn; i++)
0115   {
0116     if (ee[i] > emax)
0117     {
0118       emax = ee[i];
0119       iy0 = iyy[i];
0120       iz0 = izz[i];
0121     }
0122   }
0123
0124   if (emax <= 0) return -1;
0125
0126   int id;
0127   float etot = 0;
0128   float sz = 0;
0129   float sy = 0;
0130
0131   for (int idz = -DXY; idz <= DXY; idz++)
0132   {
0133     for (int idy = -DXY; idy <= DXY; idy++)
0134     {
0135       id = GetTowerID(iy0 + idy, iz0 + idz, nn, iyy, izz, ee);
0136       if (id >= 0)
0137       {
0138         etot += ee[id];
0139         sz += ee[id] * (iz0 + idz);
0140         sy += ee[id] * (iy0 + idy);
0141       }
0142     }
0143   }
0144   float zcg = sz / etot;  // Here cg allowed to be out of range
0145   float ycg = sy / etot;
0146   int iz0cg = int(zcg + 0.5);
0147   int iy0cg = int(ycg + 0.5);
0148   float ddz = fabs(zcg - iz0cg);
0149   float ddy = fabs(ycg - iy0cg);
0150
0151   int isz = 1;
0152   if (zcg - iz0cg < 0) isz = -1;
0153   int isy = 1;
0154   if (ycg - iy0cg < 0) isy = -1;
0155
0156   // 4 central towers: 43
0157   //                   12
0158   // Tower 1 - central one
0159   float e1, e2, e3, e4;
0160   e1 = e2 = e3 = e4 = 0;
0161   id = GetTowerID(iy0cg, iz0cg, nn, iyy, izz, ee);
0162   if (id >= 0) e1 = ee[id];
0163   id = GetTowerID(iy0cg, iz0cg + isz, nn, iyy, izz, ee);
0164   if (id >= 0) e2 = ee[id];
0165   id = GetTowerID(iy0cg + isy, iz0cg + isz, nn, iyy, izz, ee);
0166   if (id >= 0) e3 = ee[id];
0167   id = GetTowerID(iy0cg + isy, iz0cg, nn, iyy, izz, ee);
0168   if (id >= 0) e4 = ee[id];
0169
0170   float e1t = (e1 + e2 + e3 + e4) / etot;
0171   float e2t = (e1 + e2 - e3 - e4) / etot;
0172   float e3t = (e1 - e2 - e3 + e4) / etot;
0173   float e4t = (e3) / etot;
0174   //  float e5t = (e2+e4)/etot;
0175
0176   float rr = sqrt((0.5 - ddz) * (0.5 - ddz) + (0.5 - ddy) * (0.5 - ddy));
0177
0178   float c1, c2, c11;
0179
0180   float logE = log(etot);
0181
0182   // e1 energy is the most effective for PID if properly tuned !
0183   // Discrimination power is very sensitive to paramter c1: the bigger it is
0184   // the better discrimination;
0185   c1 = 0.95;
0186   c2 = 0.0066364 * logE + 0.00466667;
0187   if (c2 < 0) c2 = 0;
0188   float e1p = c1 - c2 * rr * rr;
0189   c1 = 0.034 - 0.01523 * logE + 0.0029 * logE * logE;
0190   float err1 = c1;
0191
0192   // For e2
0193   c1 = 0.00844086 + 0.00645359 * logE - 0.00119381 * logE * logE;
0194   if (etot > 15) c1 = 0.00844086 + 0.00645359 * log(15.) - 0.00119381 * log(15.) * log(15.);  // Const at etot>15GeV
0195   if (c1 < 0) c1 = 0;
0196   c2 = 3.9;                                                      // Fixed
0197   float e2p = sqrt(c1 + 0.25 * c2) - sqrt(c1 + c2 * ddy * ddy);  // =0 at ddy=0.5
0198
0199   c1 = 0.0212333 + 0.0420473 / etot;
0200   c2 = 0.090;  // Fixed
0201   float err2 = c1 + c2 * ddy;
0202   if (ddy > 0.3) err2 = c1 + c2 * 0.3;  // Const at ddy>0.3
0203
0204   // For e3
0205   c1 = 0.0107857 + 0.0056801 * logE - 0.000892016 * logE * logE;
0206   if (etot > 15) c1 = 0.0107857 + 0.0056801 * log(15.) - 0.000892016 * log(15.) * log(15.);  // Const at etot>15GeV
0207   if (c1 < 0) c1 = 0;
0208   c2 = 3.9;                                                      // Fixed
0209   float e3p = sqrt(c1 + 0.25 * c2) - sqrt(c1 + c2 * ddz * ddz);  // =0 at ddz=0.5
0210
0211   //  c1 = 0.0200 + 0.042/etot;
0212   c1 = 0.0167 + 0.058 / etot;
0213   c2 = 0.090;  // Fixed
0214   float err3 = c1 + c2 * ddz;
0215   if (ddz > 0.3) err3 = c1 + c2 * 0.3;  // Const at ddz>0.3
0216
0217   // For e4
0218   float e4p = 0.25 - 0.668 * rr + 0.460 * rr * rr;
0219   c11 = 0.171958 + 0.0142421 * logE - 0.00214827 * logE * logE;
0220   //  c11 = 0.171085 + 0.0156215*logE - -0.0025809*logE*logE;
0221   float err4 = 0.102 - 1.43 * c11 * rr + c11 * rr * rr;  // Min is set to x=1.43/2.
0222   err4 *= 1.1;
0223
0224   chi2 = 0.;
0225   chi2 += (e1p - e1t) * (e1p - e1t) / err1 / err1;
0226   chi2 += (e2p - e2t) * (e2p - e2t) / err2 / err2;
0227   chi2 += (e3p - e3t) * (e3p - e3t) / err3 / err3;
0228   chi2 += (e4p - e4t) * (e4p - e4t) / err4 / err4;
0229   ndf = 4;
0230
0231   //  chi2 /= 1.1;
0232   float prob = TMath::Prob(chi2, ndf);
0233
0234   return prob;
0235 }
0236 */
0237
0238 void BEmcRecCEMC::CorrectShowerDepth(int ix, int iy, float E, float xA, float yA, float zA, float& xC, float& yC, float& zC)
0239 {
0240   if (!m_UseCorrectShowerDepth)
0241   {
0242     xC = xA;
0243     yC = yA;
0244     zC = zA;
0245     return;
0246   }
0247
0248   float logE = log(0.1);
0249   if (E > 0.1)
0250   {
0251     logE = std::log(E);
0252   }
0253
0254   float phi = 0;
0255   // Rotate by phi (towers are tilted by a fixed angle in phi by ~9 deg?)
0256   // Just tuned from sim data
0257   if (m_UseDetailedGeometry)
0258   {
0259     phi = 0.0033 - 0.0010 * logE;
0260   }
0261   else
0262   {
0263     phi = 0.0016 - 0.0010 * logE;
0264   }
0265   xC = xA * std::cos(phi) - yA * std::sin(phi);
0266   yC = xA * std::sin(phi) + yA * std::cos(phi);
0267
0268   // Correction in z
0269   float rA = std::sqrt((xA * xA) + (yA * yA));
0270   float theta_twr;
0271   if (m_UseDetailedGeometry)
0272   {
0273     // Read the angle right from the detailed RawTowerGeom objects
0274     TowerGeom geom0;
0275     GetTowerGeometry(ix, iy, geom0);
0276     theta_twr = M_PI / 2 + geom0.rotX;
0277   }
0278   else
0279   {
0280     // Use the approximate default tower angles array.
0281     theta_twr = M_PI / 2 - angles[iy];
0282   }
0283
0284   float theta_tr = std::atan2(zA - fVz, rA);
0285
0286   // Shower CG in long. direction
0287   // Different tuning for the approximate and detailed geometry
0288   float L = 0;
0289   if (m_UseDetailedGeometry)
0290   {
0291     L = -2.67787 + 0.924138 * logE;
0292   }
0293   else
0294   {
0295     L = -1.79968 + 0.837322 * logE;
0296   }
0297   float dz = L * std::sin(theta_tr - theta_twr) / std::cos(theta_twr);
0298
0299   if (!m_UseDetailedGeometry)
0300   {
0301     // Not strictly speaking a "shower depth correction" but rather a baseline correction
0302     // The factor 0.10 accounts for the fact that the approximate geometry
0303     // is projected at 93.5cm, which is roughly 90 % of the actual average EMCal radius (~ 105 cm)
0304     dz -= fVz * 0.10;
0305   }
0306
0307   zC = zA - dz;
0308
0309   // zvtx-dependent pseudorapidity-correction
0310   // At large zvtx (> 20 cm), the sawtooth shape of the EMCal is no longer projective wrt collision point,
0311   // it introduces a bias which can be approximately corrected with the linear formula below
0312   if (m_UseDetailedGeometry && std::abs(fVz) > 20)
0313   {
0314     float eta = std::asinh((zC - fVz) / rA);
0315     if (fVz > 0)
0316     {
0317       eta -= slopes_[iy] * (fVz - 20);
0318     }
0319     else
0320     {
0321       eta -= slopes_[fNy - iy - 1] * (fVz + 20);
0322     }
0323     zC = fVz + rA * std::sinh(eta);
0324   }
0325
0326   return;
0327 }
0328
0329 void BEmcRecCEMC::CorrectEnergy(float Energy, float /*x*/, float /*y*/,
0330                                 float& Ecorr)
0331 {
0332   // Corrects the EM Shower Energy for attenuation in fibers and
0333   // long energy leakage
0334   //
0335   // (x,y) - impact position (cm) in Sector frame
0336   /*
0337   float sinT;
0338   float att, leak, corr;
0339   const float leakPar = 0.0033; // parameter from fit
0340   const float attPar = 120; // Attenuation in module (cm)
0341   const float X0 = 2; // radiation length (cm)
0342
0343   *Ecorr = Energy;
0344   if( Energy < 0.01 ) return;
0345
0346   GetImpactAngle(x, y, &sinT); // sinT not used so far
0347   leak = 2-sqrt(1+leakPar*log(1+Energy)*log(1+Energy));
0348   att = exp(log(Energy)*X0/attPar);
0349   corr = leak*att;
0350   *Ecorr = Energy/corr;
0351   */
0352   Ecorr = Energy;
0353 }
0354
0355 void BEmcRecCEMC::CorrectECore(float Ecore, float /*x*/, float /*y*/, float& Ecorr)
0356 {
0357   // Corrects the EM Shower Core Energy for attenuation in fibers,
0358   // long energy leakage and angle dependance
0359   //
0360   // (x,y) - impact position (cm) in Sector frame
0361
0362   const float c0 = 0.950;  // For no threshold
0363   Ecorr = Ecore / c0;
0364 }
0365
0366 void BEmcRecCEMC::CorrectPosition(float Energy, float x, float y,
0367                                   float& xc, float& yc)
0368 {
0369   // Corrects the Shower Center of Gravity for the systematic shift due to
0370   // limited tower size
0371   //
0372   // Everything here is in tower units.
0373   // (x,y) - CG position, (xc,yc) - corrected position
0374
0375   float bx;
0376   float by;
0377   float x0;
0378   float y0;
0379   int ix0;
0380   int iy0;
0381
0382   if (!m_UseCorrectPosition)
0383   {
0384     xc = x;
0385     yc = y;
0386     return;
0387   }
0388
0389   if (Energy < 0.01)
0390   {
0391     return;
0392   }
0393
0394   bx = 0.15;
0395   by = 0.15;
0396
0397   x0 = x;
0398   ix0 = EmcCluster::lowint(x0 + 0.5);
0399
0400   if (EmcCluster::ABS(x0 - ix0) <= 0.5)
0401   {
0402     x0 = ix0 + bx * asinh(2. * (x0 - ix0) * sinh(0.5 / bx));
0403   }
0404   else
0405   {
0406     x0 = x;
0407     std::cout << "????? Something wrong in BEmcRecCEMC::CorrectPosition: x = "
0408               << x << " dx = " << x0 - ix0 << std::endl;
0409   }
0410
0411   // Correct for phi bias within module of 8 towers
0412 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
0413   int ix8 = int(x + 0.5) / 8;
0414   float x8 = x + 0.5 - (ix8 * 8) - 4;  // from -4 to +4
0415   float dx = 0;
0416   if (m_UseDetailedGeometry)
0417   {
0418     // Don't know why there is a different factor for each tower of the sector
0419     // Just tuned from MC
0420     int local_ix8 = int(x+0.5) - ix8 * 8;
0421     dx = factor_[local_ix8] * x8 / 4.;
0422   }
0423   else
0424   {
0425     dx = 0.10 * x8 / 4.;
0426     if (std::fabs(x8) > 3.3)
0427     {
0428       dx = 0;  // Don't correct near the module edge
0429     }
0430   }
0431
0432   xc = x0 - dx;
0433   while (xc < -0.5)
0434   {
0435     xc += float(fNx);
0436   }
0437   while (xc >= fNx - 0.5)
0438   {
0439     xc -= float(fNx);
0440   }
0441
0442   y0 = y;
0443   iy0 = EmcCluster::lowint(y0 + 0.5);
0444
0445   if (EmcCluster::ABS(y0 - iy0) <= 0.5)
0446   {
0447     y0 = iy0 + by * asinh(2. * (y0 - iy0) * sinh(0.5 / by));
0448   }
0449   else
0450   {
0451     y0 = y;
0452     std::cout << "????? Something wrong in BEmcRecCEMC::CorrectPosition: y = "
0453               << y << "dy = " << y0 - iy0 << std::endl;
0454   }
0455   yc = y0;
0456 }