ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BEmcRecCEMC.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BEmcRecCEMC.cc
1 #include "BEmcRecCEMC.h"
2 
3 #include "BEmcCluster.h"
4 #include "BEmcProfile.h"
5 
6 #include <cmath>
7 #include <iostream>
8 
9 using namespace std;
10 
12 // : _emcprof(nullptr)
13 {
14  Name("BEmcRecCEMC");
15  SetCylindricalGeometry();
16 }
17 
19 {
20  // you can delete null pointers
21  // delete _emcprof;
22 }
23 
24 void BEmcRecCEMC::LoadProfile(const string& fname)
25 {
26  // cout << "Infor from BEmcRecCEMC::LoadProfile(): no external file used for shower profile evaluation in CEMC" << endl;
27  _emcprof = new BEmcProfile(fname);
28 }
29 
30 
31 void BEmcRecCEMC::GetImpactThetaPhi(float xg, float yg, float zg, float& theta, float& phi)
32 {
33  theta = 0;
34  phi = 0;
35 
36  // float theta = atan(sqrt(xg*xg + yg*yg)/fabs(zg-fVz));
37  float rg = sqrt(xg*xg+yg*yg);
38  float theta_twr;
39  if( fabs(zg)<=15 ) theta_twr = 0;
40  else if( zg>15 ) theta_twr = atan2(zg-15,rg);
41  else theta_twr = atan2(zg+15,rg);
42  float theta_tr = atan2(zg-fVz,rg);
43  theta = fabs(theta_tr - theta_twr);
44  // phi = atan2(yg,xg);
45 }
46 
47 /*
48 float BEmcRecCEMC::GetProb(vector<EmcModule> HitList, float ecl, float xg, float yg, float zg, float& chi2, int& ndf)
49 {
50  chi2 = 0;
51  ndf = 0;
52  float prob = -1;
53 
54  // float theta = atan(sqrt(xg*xg + yg*yg)/fabs(zg-fVz));
55  float rg = sqrt(xg * xg + yg * yg);
56  float theta_twr;
57  if (fabs(zg) <= 15)
58  theta_twr = 0;
59  else if (zg > 15)
60  theta_twr = atan2(zg - 15, rg);
61  else
62  theta_twr = atan2(zg + 15, rg);
63  float theta_tr = atan2(zg - fVz, rg);
64  float theta = fabs(theta_tr - theta_twr);
65 
66  float phi = atan2(yg, xg);
67  if (_emcprof != nullptr) prob = _emcprof->GetProb(&HitList, fNx, ecl, theta, phi);
68 
69  return prob;
70 }
71 */
72 /*
73 float BEmcRecCEMC::GetProb(vector<EmcModule> HitList, float et, float xg, float yg, float zg, float& chi2, int& ndf)
74 // et, xg, yg, zg not used here
75 {
76  const float thresh = 0.01;
77  const int DXY = 3; // 2 is for 5x5 matrix; 3 for 7x7 matrix
78  const int Nmax = 1000;
79  float ee[Nmax];
80  int iyy[Nmax];
81  int izz[Nmax];
82 
83  int ich;
84  vector<EmcModule>::iterator ph = HitList.begin();
85 
86  chi2 = 0;
87  ndf = 0;
88 
89  int nn = 0;
90 
91  while (ph != HitList.end())
92  {
93  ee[nn] = ph->amp;
94  if (ee[nn] > thresh)
95  {
96  ich = ph->ich;
97  izz[nn] = ich % fNx;
98  iyy[nn] = ich / fNx;
99  nn++;
100  if (nn >= Nmax)
101  {
102  cout << "BEmcRec::GetProb: Cluster size is too big. Skipping the rest of the towers" << endl;
103  break;
104  }
105  } // if( ee[nn]
106  ++ph;
107  } // while( ph
108 
109  if (nn <= 0) return -1;
110 
111  int iy0 = -1, iz0 = -1;
112  float emax = 0;
113 
114  for (int i = 0; i < nn; i++)
115  {
116  if (ee[i] > emax)
117  {
118  emax = ee[i];
119  iy0 = iyy[i];
120  iz0 = izz[i];
121  }
122  }
123 
124  if (emax <= 0) return -1;
125 
126  int id;
127  float etot = 0;
128  float sz = 0;
129  float sy = 0;
130 
131  for (int idz = -DXY; idz <= DXY; idz++)
132  {
133  for (int idy = -DXY; idy <= DXY; idy++)
134  {
135  id = GetTowerID(iy0 + idy, iz0 + idz, nn, iyy, izz, ee);
136  if (id >= 0)
137  {
138  etot += ee[id];
139  sz += ee[id] * (iz0 + idz);
140  sy += ee[id] * (iy0 + idy);
141  }
142  }
143  }
144  float zcg = sz / etot; // Here cg allowed to be out of range
145  float ycg = sy / etot;
146  int iz0cg = int(zcg + 0.5);
147  int iy0cg = int(ycg + 0.5);
148  float ddz = fabs(zcg - iz0cg);
149  float ddy = fabs(ycg - iy0cg);
150 
151  int isz = 1;
152  if (zcg - iz0cg < 0) isz = -1;
153  int isy = 1;
154  if (ycg - iy0cg < 0) isy = -1;
155 
156  // 4 central towers: 43
157  // 12
158  // Tower 1 - central one
159  float e1, e2, e3, e4;
160  e1 = e2 = e3 = e4 = 0;
161  id = GetTowerID(iy0cg, iz0cg, nn, iyy, izz, ee);
162  if (id >= 0) e1 = ee[id];
163  id = GetTowerID(iy0cg, iz0cg + isz, nn, iyy, izz, ee);
164  if (id >= 0) e2 = ee[id];
165  id = GetTowerID(iy0cg + isy, iz0cg + isz, nn, iyy, izz, ee);
166  if (id >= 0) e3 = ee[id];
167  id = GetTowerID(iy0cg + isy, iz0cg, nn, iyy, izz, ee);
168  if (id >= 0) e4 = ee[id];
169 
170  float e1t = (e1 + e2 + e3 + e4) / etot;
171  float e2t = (e1 + e2 - e3 - e4) / etot;
172  float e3t = (e1 - e2 - e3 + e4) / etot;
173  float e4t = (e3) / etot;
174  // float e5t = (e2+e4)/etot;
175 
176  float rr = sqrt((0.5 - ddz) * (0.5 - ddz) + (0.5 - ddy) * (0.5 - ddy));
177 
178  float c1, c2, c11;
179 
180  float logE = log(etot);
181 
182  // e1 energy is the most effective for PID if properly tuned !
183  // Discrimination power is very sensitive to paramter c1: the bigger it is
184  // the better discrimination;
185  c1 = 0.95;
186  c2 = 0.0066364 * logE + 0.00466667;
187  if (c2 < 0) c2 = 0;
188  float e1p = c1 - c2 * rr * rr;
189  c1 = 0.034 - 0.01523 * logE + 0.0029 * logE * logE;
190  float err1 = c1;
191 
192  // For e2
193  c1 = 0.00844086 + 0.00645359 * logE - 0.00119381 * logE * logE;
194  if (etot > 15) c1 = 0.00844086 + 0.00645359 * log(15.) - 0.00119381 * log(15.) * log(15.); // Const at etot>15GeV
195  if (c1 < 0) c1 = 0;
196  c2 = 3.9; // Fixed
197  float e2p = sqrt(c1 + 0.25 * c2) - sqrt(c1 + c2 * ddy * ddy); // =0 at ddy=0.5
198 
199  c1 = 0.0212333 + 0.0420473 / etot;
200  c2 = 0.090; // Fixed
201  float err2 = c1 + c2 * ddy;
202  if (ddy > 0.3) err2 = c1 + c2 * 0.3; // Const at ddy>0.3
203 
204  // For e3
205  c1 = 0.0107857 + 0.0056801 * logE - 0.000892016 * logE * logE;
206  if (etot > 15) c1 = 0.0107857 + 0.0056801 * log(15.) - 0.000892016 * log(15.) * log(15.); // Const at etot>15GeV
207  if (c1 < 0) c1 = 0;
208  c2 = 3.9; // Fixed
209  float e3p = sqrt(c1 + 0.25 * c2) - sqrt(c1 + c2 * ddz * ddz); // =0 at ddz=0.5
210 
211  // c1 = 0.0200 + 0.042/etot;
212  c1 = 0.0167 + 0.058 / etot;
213  c2 = 0.090; // Fixed
214  float err3 = c1 + c2 * ddz;
215  if (ddz > 0.3) err3 = c1 + c2 * 0.3; // Const at ddz>0.3
216 
217  // For e4
218  float e4p = 0.25 - 0.668 * rr + 0.460 * rr * rr;
219  c11 = 0.171958 + 0.0142421 * logE - 0.00214827 * logE * logE;
220  // c11 = 0.171085 + 0.0156215*logE - -0.0025809*logE*logE;
221  float err4 = 0.102 - 1.43 * c11 * rr + c11 * rr * rr; // Min is set to x=1.43/2.
222  err4 *= 1.1;
223 
224  chi2 = 0.;
225  chi2 += (e1p - e1t) * (e1p - e1t) / err1 / err1;
226  chi2 += (e2p - e2t) * (e2p - e2t) / err2 / err2;
227  chi2 += (e3p - e3t) * (e3p - e3t) / err3 / err3;
228  chi2 += (e4p - e4t) * (e4p - e4t) / err4 / err4;
229  ndf = 4;
230 
231  // chi2 /= 1.1;
232  float prob = TMath::Prob(chi2, ndf);
233 
234  return prob;
235 }
236 */
237 
238 void BEmcRecCEMC::CorrectShowerDepth(float E, float xA, float yA, float zA, float& xC, float& yC, float& zC)
239 {
240  /*
241  xC = xA;
242  yC = yA;
243  zC = zA;
244  */
245 
246  float logE = log(0.1);
247  if( E>0.1 ) logE = log(E);
248 
249  // Rotate by phi (towers are tilted by a fixed angle in phi by ~9 deg?)
250  // Just tuned from sim data
251  float phi = 0.002 - 0.001*logE;
252  xC = xA*cos(phi) - yA*sin(phi);
253  yC = xA*sin(phi) + yA*cos(phi);
254 
255  // Correction in z
256  // Just tuned for sim data ... don't fully understand why it works like that
257  float rA = sqrt(xA*xA + yA*yA);
258  // float theta_twr = GetTowerTheta(xA,yA,zA);
259  float theta_twr;
260  if( fabs(zA)<=15 ) theta_twr = 0;
261  else if( zA>15 ) theta_twr = atan2(zA-15,rA);
262  else theta_twr = atan2(zA+15,rA);
263 
264  float theta_tr = atan2(zA-fVz,rA);
265  float L = -1.3 + 0.7*logE; // Shower CG in long. direction
266  float dz = L*sin(theta_tr-theta_twr)/cos(theta_twr);
267 
268  dz -= fVz*0.10;
269 
270  zC = zA - dz;
271 
272  return;
273 }
274 
275 void BEmcRecCEMC::CorrectEnergy(float Energy, float /*x*/, float /*y*/,
276  float& Ecorr)
277 {
278  // Corrects the EM Shower Energy for attenuation in fibers and
279  // long energy leakage
280  //
281  // (x,y) - impact position (cm) in Sector frame
282  /*
283  float sinT;
284  float att, leak, corr;
285  const float leakPar = 0.0033; // parameter from fit
286  const float attPar = 120; // Attenuation in module (cm)
287  const float X0 = 2; // radiation length (cm)
288 
289  *Ecorr = Energy;
290  if( Energy < 0.01 ) return;
291 
292  GetImpactAngle(x, y, &sinT); // sinT not used so far
293  leak = 2-sqrt(1+leakPar*log(1+Energy)*log(1+Energy));
294  att = exp(log(Energy)*X0/attPar);
295  corr = leak*att;
296  *Ecorr = Energy/corr;
297  */
298  Ecorr = Energy;
299 }
300 
301 void BEmcRecCEMC::CorrectECore(float Ecore, float /*x*/, float /*y*/, float& Ecorr)
302 {
303  // Corrects the EM Shower Core Energy for attenuation in fibers,
304  // long energy leakage and angle dependance
305  //
306  // (x,y) - impact position (cm) in Sector frame
307 
308  const float c0 = 0.950; // For no threshold
309  Ecorr = Ecore / c0;
310 }
311 
312 void BEmcRecCEMC::CorrectPosition(float Energy, float x, float y,
313  float& xc, float& yc)
314 {
315  // Corrects the Shower Center of Gravity for the systematic shift due to
316  // limited tower size
317  //
318  // Everything here is in tower units.
319  // (x,y) - CG position, (xc,yc) - corrected position
320 
321  float xZero, yZero, bx, by;
322  float t, x0, y0;
323  int ix0, iy0;
324 
325  xc = x;
326  yc = y;
327 
328  if (Energy < 0.01) return;
329  /*
330  float xA, yA, zA;
331  Tower2Global(Energy, x, y, xA, yA, zA);
332  zA -= fVz;
333  float sinTx = xA / sqrt(xA * xA + zA * zA);
334  float sinTy = yA / sqrt(yA * yA + zA * zA);
335  */
336  float sinTx = 0;
337  float sinTy = 0;
338 
339  float sin2Tx = sinTx * sinTx;
340  float sin2Ty = sinTy * sinTy;
341 
342  if (sinTx > 0)
343  xZero = -0.417 * sinTx - 1.500 * sin2Tx;
344  else
345  xZero = -0.417 * sinTx + 1.500 * sin2Tx;
346 
347  if (sinTy > 0)
348  yZero = -0.417 * sinTy - 1.500 * sin2Ty;
349  else
350  yZero = -0.417 * sinTy + 1.500 * sin2Ty;
351 
352  t = 0.98 + 0.98 * sqrt(Energy);
353  bx = 0.15 + t * sin2Tx;
354  by = 0.15 + t * sin2Ty;
355 
356  x0 = x + xZero;
357  ix0 = EmcCluster::lowint(x0 + 0.5);
358 
359  if (EmcCluster::ABS(x0 - ix0) <= 0.5)
360  {
361  x0 = (ix0 - xZero) + bx * asinh(2. * (x0 - ix0) * sinh(0.5 / bx));
362  }
363  else
364  {
365  x0 = x;
366  cout << "????? Something wrong in BEmcRecCEMC::CorrectPosition: x = "
367  << x << " dx = " << x0 - ix0 << endl;
368  }
369 
370  // Correct for phi bias within module of 8 towers
371  int ix8 = int(x + 0.5) / 8;
372  float x8 = x + 0.5 - ix8 * 8 - 4; // from -4 to +4
373  float dx = 0.10 * x8 / 4.;
374  if (fabs(x8) > 3.3) dx = 0; // Don't correct near the module edge
375  // dx = 0;
376 
377  xc = x0 - dx;
378  while (xc < -0.5) xc += float(fNx);
379  while (xc >= fNx - 0.5) xc -= float(fNx);
380 
381  y0 = y + yZero;
382  iy0 = EmcCluster::lowint(y0 + 0.5);
383 
384  if (EmcCluster::ABS(y0 - iy0) <= 0.5)
385  {
386  y0 = (iy0 - yZero) + by * asinh(2. * (y0 - iy0) * sinh(0.5 / by));
387  }
388  else
389  {
390  y0 = y;
391  cout << "????? Something wrong in BEmcRecCEMC::CorrectPosition: y = "
392  << y << "dy = " << y0 - iy0 << endl;
393  }
394  yc = y0;
395 }
396 
397 /*
398 void BEmcRecCEMC::CorrectPosition(float Energy, float x, float y,
399  float* pxc, float* pyc)
400 {
401  // Corrects the Shower Center of Gravity for the systematic error due to
402  // the limited tower size and angle shift
403  //
404  // Everything here is in cell units.
405  // (x,y) - CG position, (*pxc,*pyc) - corrected position
406 
407  float xShift, yShift, xZero, yZero, bx, by;
408  float t, x0, y0;
409  int ix0, iy0;
410  // int signx, signy;
411 
412  const float Xrad = 0.3; // !!!!! Need to put correct value
413  const float Remc = 90.; // EMCal inner radius. !!!!! Should be obtained from geometry container
414 
415  *pxc = x;
416  *pyc = y;
417  // return;
418 
419  SetProfileParameters(0, Energy, x, y);
420  // if( fSinTx >= 0 ) signx = 1;
421  // else signx = -1;
422  // if( fSinTy >= 0 ) signy = 1;
423  // else signy = -1;
424  t = 5.0 + 1.0 * log(Energy); // In Rad Length units
425  t *= (Xrad / Remc / GetModSizex()); // !!!!!
426  xShift = t * fSinTx;
427  yShift = t * fSinTy;
428  // xZero=xShift-(0.417*EmcCluster::ABS(fSinTx)+1.500*fSinTx*fSinTx)*signx;
429  // yZero=yShift-(0.417*EmcCluster::ABS(fSinTy)+1.500*fSinTy*fSinTy)*signy;
430  xZero = xShift; // ...Somehow this works better !!!!!
431  yZero = yShift; // ...Somehow this works better !!!!!
432  t = 0.98 + 0.98 * sqrt(Energy); // !!!!! Still from PHENIX
433  bx = 0.15 + t * fSinTx * fSinTx;
434  by = 0.15 + t * fSinTy * fSinTy;
435 
436  x0 = x;
437  x0 = x0 - xShift + xZero;
438  ix0 = EmcCluster::lowint(x0 + 0.5);
439  if (EmcCluster::ABS(x0 - ix0) <= 0.5)
440  {
441  x0 = (ix0 - xZero) + bx * asinh(2. * (x0 - ix0) * sinh(0.5 / bx));
442  *pxc = x0;
443  }
444  else
445  {
446  *pxc = x - xShift;
447  cout << "????? Something wrong in CorrectPosition: x = "
448  << x << " dx = " << x0 - ix0 << endl;
449  }
450 
451  y0 = y;
452  y0 = y0 - yShift + yZero;
453  iy0 = EmcCluster::lowint(y0 + 0.5);
454  if (EmcCluster::ABS(y0 - iy0) <= 0.5)
455  {
456  y0 = (iy0 - yZero) + by * asinh(2. * (y0 - iy0) * sinh(0.5 / by));
457  *pyc = y0;
458  }
459  else
460  {
461  *pyc = y - yShift;
462  cout << "????? Something wrong in CorrectPosition: y = "
463  << y << " dy = " << y0 - iy << endl;
464  }
465 }
466 */