ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BEmcRecEEMC.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BEmcRecEEMC.cc
1 #include "BEmcRecEEMC.h"
2 
3 #include "BEmcCluster.h"
4 #include "BEmcProfile.h"
5 
6 #include <calobase/RawTowerDefs.h>
7 
8 #include <cmath>
9 #include <iostream>
10 
12 {
13  Name("BEmcRecEEMC");
15 }
16 
17 void BEmcRecEEMC::LoadProfile(const std::string& fname)
18 {
19  // std::cout << "Info from BEmcRecEEMC::LoadProfile(): no shower profile evaluation is defined yet for EEMC" << std::endl;
20  _emcprof = new BEmcProfile(fname);
21 }
22 
23 void BEmcRecEEMC::GetImpactThetaPhi(float xg, float yg, float zg, float& theta, float& phi)
24 {
25  theta = atan(sqrt(xg*xg + yg*yg)/fabs(zg-fVz));
26  phi = atan2(yg,xg);
27 }
28 
29 /*
30 float BEmcRecEEMC::GetProb(vector<EmcModule> HitList, float ecl, float xg, float yg, float zg, float& chi2, int& ndf)
31 // ecl, xg, yg, zg not used here
32 {
33  chi2 = 0;
34  ndf = 0;
35  float prob = -1;
36 
37  float theta = atan(sqrt(xg*xg + yg*yg)/fabs(zg-fVz));
38  float phi = atan2(yg,xg);
39  if( _emcprof != nullptr ) prob = _emcprof->GetProb(&HitList,fNx,ecl,theta,phi);
40 
41  return prob;
42 }
43 */
44 
45 void BEmcRecEEMC::CorrectShowerDepth(float E, float xA, float yA, float zA, float& xC, float& yC, float& zC)
46 {
47  // DZ, D and X0 should be negative for -Z direction
48  // Crystal:[-10., -6.1, -0.91], Sci_glass:[-20., -16.76, -2.5]
49  // Original one:[-9, -6.1, -0.91] also for the lead tungsten
50 
51  int C_ID = BEmcRec::GetCalotype();
52  // std::cout << "Success pass data(ID): " << C_ID << std::endl;
53 
54 
55  if( (C_ID == RawTowerDefs::EEMC) || (C_ID == RawTowerDefs::EEMC_crystal) )
56  {
57  const float DZ = -0.5 * BEmcRec::GetScinSize(); // in cm, tower half length
58  const float D = -6.1; // in cm, shower depth at 1 GeV relative to tower face; obtained from GEANT
59  const float X0 = -0.91; // in cm; obtained from GEANT (should be ~ rad length)
60 
61  // std::cout << "Success pass data(size): " << DZ << std::endl << std::endl;
62 
63  float logE = log(0.1);
64  if (E > 0.1) logE = log(E);
65  float zV = zA - fVz;
66  float cosT = fabs(zV) / sqrt(xA * xA + yA * yA + zV * zV);
67 
68  zC = (zA - DZ) + (D + X0 * logE) * cosT; //Only the shower depth corrected
69  // zC = zA; // !!!!!
70 
71  xC = xA; // Keep the x and y the same. The x and y correction is in another code
72  yC = yA;
73  }
74  else if( C_ID == RawTowerDefs::EEMC_glass )
75  {
76  const float DZ = -0.5 * BEmcRec::GetScinSize(); // in cm, tower half length
77  const float D = -15.25; // in cm, shower depth at 1 GeV relative to tower face; obtained from GEANT
78  const float X0 = -2.275; // in cm; obtained from GEANT (should be ~ rad length)
79 
80  // std::cout << "Success pass data(size): " << DZ << std::endl << std::endl;
81 
82  float logE = log(0.1);
83  if (E > 0.1) logE = log(E);
84  float zV = zA - fVz;
85  float cosT = fabs(zV) / sqrt(xA * xA + yA * yA + zV * zV);
86 
87  zC = (zA - DZ) + (D + X0 * logE) * cosT; //Only the shower depth corrected
88  // zC = zA; // !!!!!
89 
90  xC = xA; // Keep the x and y the same. The x and y correction is in another code
91  yC = yA;
92  }
93 
94 }
95 
96 void BEmcRecEEMC::CorrectEnergy(float Energy, float /*x*/, float /*y*/,
97  float& Ecorr)
98 {
99  Ecorr = Energy;
100 }
101 
102 void BEmcRecEEMC::CorrectECore(float Ecore, float /*x*/, float /*y*/, float& Ecorr)
103 {
104  // Corrects the EM Shower Core Energy for attenuation in fibers,
105  // long energy leakage and angle dependance
106  //
107  // (x,y) - shower CG in tower units (not projected anywhere!)
108 
109  Ecorr = Ecore;
110 }
111 
112 float BEmcRecEEMC::GetImpactAngle(float /*e*/, float /*x*/, float /*y*/)
113 // Get impact angle, (x,y) - position in Sector frame (cm)
114 {
115  /*
116  float xVert, yVert, zVert;
117  float vx, vy, vz;
118 
119  GlobalToSector( fVx, fVy, fVz, &xVert, &yVert, &zVert );
120  vz = -zVert;
121  vy = y - yVert;
122  vx = x - xVert;
123  // From this point X coord in sector frame is Z coord in Global Coord System !!!
124  *sinT = sqrt((vx*vx+vy*vy)/(vx*vx+vy*vy+vz*vz));
125  */
126  return 0;
127 }
128 
129 void BEmcRecEEMC::CorrectPosition(float Energy, float x, float y,
130  float& xc, float& yc)
131 {
132  // Corrects the Shower Center of Gravity for the systematic shift due to
133  // the limited tower size
134  //
135  // Everything here is in tower units.
136  // (x,y) - CG position, (xc,yc) - corrected position
137 
138  float xZero, yZero, bx, by;
139  float t, x0, y0;
140  int ix0, iy0;
141 
142  xc = x;
143  yc = y;
144  // return;
145 
146  if (Energy < 0.01) return;
147 
148  float xA, yA, zA;
149  Tower2Global(Energy, x, y, xA, yA, zA);
150  zA -= fVz;
151  // float sinTx = xA / sqrt(xA * xA + zA * zA);
152  // float sinTy = yA / sqrt(yA * yA + zA * zA);
153  float sinTy = xA / sqrt(xA * xA + zA * zA); // x is second index in here
154  float sinTx = yA / sqrt(yA * yA + zA * zA);
155  float sin2Tx = sinTx * sinTx;
156  float sin2Ty = sinTy * sinTy;
157 
158  if (sinTx > 0)
159  xZero = -0.6 * sinTx - 1.500 * sin2Tx;
160  else
161  xZero = -0.6 * sinTx + 1.500 * sin2Tx;
162 
163  if (sinTy > 0)
164  yZero = -0.6 * sinTy - 1.5 * sin2Ty;
165  else
166  yZero = -0.6 * sinTy + 1.5 * sin2Ty;
167 
168  yZero = -yZero; // Because tower index in X decreases with increasing X in EEMC (y is actually x here!)
169 
170  t = 0.98 + 0.98 * sqrt(Energy);
171  t *= 0.7; // Temp correction
172  bx = 0.20 - 0.009 * log(Energy) + t * sin2Tx;
173  by = 0.20 - 0.009 * log(Energy) + t * sin2Ty;
174 
175  x0 = x + xZero;
176  ix0 = EmcCluster::lowint(x0 + 0.5);
177 
178 
179 
180  if (EmcCluster::ABS(x0 - ix0) <= 0.5)
181  {
182  x0 = (ix0 - xZero) + bx * asinh(2. * (x0 - ix0) * sinh(0.5 / bx));
183  xc = x0;
184  }
185  else
186  {
187  xc = x;
188  std::cout << "????? Something wrong in BEmcRecEEMC::CorrectPosition: x = "
189  << x << ", dx = " << x0 - ix0 << std::endl;
190  }
191 
192  y0 = y + yZero;
193  iy0 = EmcCluster::lowint(y0 + 0.5);
194 
195  if (EmcCluster::ABS(y0 - iy0) <= 0.5)
196  {
197  y0 = (iy0 - yZero) + by * asinh(2. * (y0 - iy0) * sinh(0.5 / by));
198  yc = y0;
199  }
200  else
201  {
202  yc = y;
203  std::cout << "????? Something wrong in BEmcRecEEMC::CorrectPosition: y = "
204  << y << ", dy = " << y0 - iy0 << std::endl;
205  }
206 
207 }