ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BEmcRec.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BEmcRec.cc
1 // Name: BEmcRec.h
2 // Author: A. Bazilevsky, Apr 2012
3 // Modified from EmcSectorRec.cxx and EmcScSectorRec.cxx
4 
5 // BEmcRec -- base class for sPHENIX EMCal
6 
7 // ///////////////////////////////////////////////////////////////////////////
8 
9 #include "BEmcRec.h"
10 #include "BEmcCluster.h"
11 #include "BEmcProfile.h"
12 
13 #include <TMath.h>
14 
15 #include <cstdlib>
16 #include <fstream>
17 #include <iostream>
18 #include <utility>
19 
20 // ///////////////////////////////////////////////////////////////////////////
21 // BEmcRec member functions
22 
24 {
25  fTowerGeom.clear();
26  fModules = new std::vector<EmcModule>;
27  fClusters = new std::vector<EmcCluster>;
28 }
29 
30 // ///////////////////////////////////////////////////////////////////////////
31 
33 {
34  if (fModules)
35  {
36  fModules->clear();
37  delete fModules;
38  }
39 
40  if (fClusters)
41  {
42  fClusters->clear();
43  delete fClusters;
44  }
45 
46  delete _emcprof;
47 }
48 
49 // ///////////////////////////////////////////////////////////////////////////
50 
51 void BEmcRec::LoadProfile(const std::string& /*fname*/)
52 {
53  std::cout << "Warning from BEmcRec::LoadProfile(): No acton defined for shower profile evaluation; should be defined in a detector specific module " << Name() << std::endl;
54 }
55 
56 void BEmcRec::PrintTowerGeometry(const std::string& fname)
57 {
58  std::ofstream outfile(fname);
59  if (!outfile.is_open())
60  {
61  std::cout << "Error in BEmcRec::PrintTowerGeometry(): Failed to open file "
62  << fname << std::endl;
63  return;
64  }
65  outfile << "Number of bins:" << std::endl;
66  outfile << fNx << " " << fNy << std::endl;
67  outfile << "ix iy x y z dx0 dy0 dz0 dx1 dy1 dz1" << std::endl;
68  int ich;
69  TowerGeom geom;
70  std::map<int, TowerGeom>::iterator it;
71  for (int iy = 0; iy < fNy; iy++)
72  {
73  for (int ix = 0; ix < fNx; ix++)
74  {
75  ich = iy * fNx + ix;
76  it = fTowerGeom.find(ich);
77  if (it != fTowerGeom.end())
78  {
79  geom = it->second;
80  outfile << ix << " " << iy << " " << geom.Xcenter << " "
81  << geom.Ycenter << " " << geom.Zcenter << " " << geom.dX[0] << " "
82  << geom.dY[0] << " " << geom.dZ[0] << " " << geom.dX[1] << " "
83  << geom.dY[1] << " " << geom.dZ[1] << std::endl;
84  // std::cout << "Z0: " << geom.dZ[0] << " || Z1: " << geom.dZ[1] << std::endl;
85  }
86  }
87  }
88 }
89 
90 bool BEmcRec::GetTowerGeometry(int ix, int iy, TowerGeom& geom)
91 {
92  if (ix < 0 || ix >= fNx || iy < 0 || iy >= fNy) return false;
93 
94  int ich = iy * fNx + ix;
95  std::map<int, TowerGeom>::iterator it = fTowerGeom.find(ich);
96  if (it == fTowerGeom.end()) return false;
97 
98  geom = it->second;
99  return true;
100 }
101 
102 bool BEmcRec::SetTowerGeometry(int ix, int iy, float xx, float yy, float zz)
103 {
104  if (ix < 0 || ix >= fNx || iy < 0 || iy >= fNy) return false;
105 
106  TowerGeom geom;
107  geom.Xcenter = xx;
108  geom.Ycenter = yy;
109  geom.Zcenter = zz;
110  geom.dX[0] = geom.dX[1] = 0; // These should be calculated by CompleteTowerGeometry()
111  geom.dY[0] = geom.dY[1] = 0;
112  geom.dZ[0] = geom.dZ[1] = 0;
113 
114  int ich = iy * fNx + ix;
115  fTowerGeom[ich] = geom;
116  return true;
117 }
118 
120 // Calculates tower front size from coordinates of tower center coordinates
121 {
122  if (fTowerGeom.empty() || fNx <= 0)
123  {
124  std::cout << "Error in BEmcRec::CalculateTowerSize(): Tower geometry not well setup (NX = "
125  << fNx << ")" << std::endl;
126  return false;
127  }
128 
129  const int nb = 8;
130  int idx[nb] = {0, 1, 0, -1, -1, 1, 1, -1};
131  int idy[nb] = {-1, 0, 1, 0, -1, -1, 1, 1};
132 
133  std::map<int, TowerGeom>::iterator it;
134 
135  for (it = fTowerGeom.begin(); it != fTowerGeom.end(); ++it)
136  {
137  int ich = it->first;
138  TowerGeom geom0 = it->second;
139  int ix = ich % fNx;
140  int iy = ich / fNx;
141 
142  TowerGeom geomx;
143  int inx = 0;
144 
145  while (inx < nb && (idx[inx] == 0 || !GetTowerGeometry(ix + idx[inx], iy + idy[inx], geomx))) inx++;
146  if (inx >= nb)
147  {
148  std::cout << "Error in BEmcRec::CompleteTowerGeometry(): Error when locating neighbour for (ix,iy)=("
149  << ix << "," << iy << ")" << std::endl;
150  return false;
151  }
152 
153  TowerGeom geomy;
154  int iny = 0;
155 
156  while (iny < nb && (idy[iny] == 0 || !GetTowerGeometry(ix + idx[iny], iy + idy[iny], geomy))) iny++;
157  if (iny >= nb)
158  {
159  std::cout << "Error in BEmcRec::CompleteTowerGeometry(): Error when locating neighbour for (ix,iy)=("
160  << ix << "," << iy << ")" << std::endl;
161  return false;
162  }
163 
164  geom0.dX[0] = (geomx.Xcenter - geom0.Xcenter) / float(idx[inx]);
165  geom0.dY[0] = (geomx.Ycenter - geom0.Ycenter) / float(idx[inx]);
166  geom0.dZ[0] = (geomx.Zcenter - geom0.Zcenter) / float(idx[inx]);
167  geom0.dX[1] = (geomy.Xcenter - geom0.Xcenter) / float(idy[iny]);
168  geom0.dY[1] = (geomy.Ycenter - geom0.Ycenter) / float(idy[iny]);
169  geom0.dZ[1] = (geomy.Zcenter - geom0.Zcenter) / float(idy[iny]);
170 
171  it->second = geom0;
172 
173  } // it = fTowerGeom.begin()
174 
175  return true;
176 }
177 
178 void BEmcRec::Tower2Global(float E, float xC, float yC,
179  float& xA, float& yA, float& zA)
180 // xC and yC are local position in tower units
181 // For CYL geometry (xC,yC) is actually (phiC,zC)
182 {
183  xA = 0;
184  yA = 0;
185  zA = 0;
186 
187  int ix = xC + 0.5; // tower #
188  if (ix < 0 || ix >= fNx)
189  {
190  std::cout << m_ThisName << " Error in BEmcRec::Tower2Global: wrong input x: " << ix << std::endl;
191  return;
192  }
193 
194  int iy = yC + 0.5; // tower #
195  if (iy < 0 || iy >= fNy)
196  {
197  std::cout << "Error in BEmcRec::Tower2Global: wrong input y: " << iy << std::endl;
198  return;
199  }
200 
201  // Get tower where the shower is positioned
202  TowerGeom geom0;
203 
204  if (!GetTowerGeometry(ix, iy, geom0))
205  {
206  // Weird case: cluster center of gravity outside the EMCal, take geometry from the neighbouring tower
207  int idx[4] = {1, 0, -1, 0};
208  int idy[4] = {0, 1, 0, -1};
209  int ii = 0;
210  while (ii < 4 && !GetTowerGeometry(ix + idx[ii], iy + idy[ii], geom0)) ii++;
211  if (ii >= 4)
212  {
213  std::cout << "Error in BEmcRec::Tower2Global: can not identify neighbour for tower ("
214  << ix << "," << iy << ")" << std::endl;
215  return;
216  }
217  float Xc = geom0.Xcenter - idx[ii] * geom0.dX[0] - idy[ii] * geom0.dX[1];
218  float Yc = geom0.Ycenter - idx[ii] * geom0.dY[0] - idy[ii] * geom0.dY[1];
219  float Zc = geom0.Zcenter - idx[ii] * geom0.dZ[0] - idy[ii] * geom0.dZ[1];
220  geom0.Xcenter = Xc;
221  geom0.Ycenter = Yc;
222  geom0.Zcenter = Zc;
223  }
224 
225  float xt = geom0.Xcenter + (xC - ix) * geom0.dX[0] + (yC - iy) * geom0.dX[1];
226  float yt = geom0.Ycenter + (xC - ix) * geom0.dY[0] + (yC - iy) * geom0.dY[1];
227  float zt = geom0.Zcenter + (xC - ix) * geom0.dZ[0] + (yC - iy) * geom0.dZ[1];
228 
229  CorrectShowerDepth(E, xt, yt, zt, xA, yA, zA);
230 
231  // rA = sqrt(xA*xA+yA*yA);
232  // phiA = atan2(yA, xA);
233 }
234 
235 // ///////////////////////////////////////////////////////////////////////////
236 
237 int BEmcRec::iTowerDist(int ix1, int ix2)
238 // Distrance in tower units
239 {
240  int idist = ix2 - ix1;
241  if (bCYL)
242  {
243  int idistr = fNx - abs(idist); // Always >0
244  if (idistr < abs(idist))
245  { // Then count in opposite direction
246  if (idist < 0)
247  idist = idistr;
248  else
249  idist = -idistr;
250  }
251  }
252  // std::cout << "Dist " << ix1 << " " << ix2 << ": " << idist << std::endl;
253  return idist;
254 }
255 
256 float BEmcRec::fTowerDist(float x1, float x2)
257 {
258  float dist = x2 - x1;
259  if (bCYL)
260  {
261  float distr = fNx - fabs(dist); // Always >0
262  if (distr < abs(dist))
263  { // Then count in opposite direction
264  if (dist < 0)
265  dist = distr;
266  else
267  dist = -distr;
268  }
269  }
270  return dist;
271 }
272 
273 // ///////////////////////////////////
274 
276 // Cluster search algorithm based on Lednev's one developed for GAMS.
277 // Returns number of clusters found
278 {
279  int nhit, nCl;
280  // int LenCl[fgMaxLen];
281  int* LenCl;
282  int next, ib, ie, iab, iae, last, LastCl, leng, ich;
283  int ia = 0;
284 
285  EmcModule* vv;
286  EmcModule *vhit, *vt;
287  EmcCluster Clt(this);
288  std::vector<EmcModule>::iterator ph;
289  std::vector<EmcModule> hl;
290 
291  (*fClusters).clear();
292  nhit = (*fModules).size();
293 
294  if (nhit <= 0) return 0;
295  if (nhit == 1)
296  {
297  Clt.ReInitialize((*fModules));
298  fClusters->push_back(Clt);
299  return 1;
300  }
301 
302  int MaxLen = fgMaxLen;
303  LenCl = new int[MaxLen];
304  ZeroVector(LenCl, MaxLen);
305 
306  vt = new EmcModule[nhit];
307  vhit = new EmcModule[nhit];
308 
309  ph = (*fModules).begin();
310  vv = vhit;
311  while (ph != (*fModules).end()) *vv++ = *ph++;
312 
313  qsort(vhit, nhit, sizeof(EmcModule), HitNCompare);
314 
315  nCl = 0;
316  next = 0;
317  for (ich = 1; ich < nhit + 1; ich++)
318  {
319  if (ich < nhit) ia = vhit[ich].ich;
320 
321  // New subcluster
322  //
323  if ((ia - vhit[ich - 1].ich > 1) // the beginning of new subcluster
324  || (ich >= nhit) // just finish defining last sub-cluster
325  || (ia - ia / fNx * fNx == 0))
326  { // new raw -> new subcluster
327 
328  ib = next;
329  ie = ich - 1;
330  next = ich;
331  if (nCl >= MaxLen)
332  {
333  // delete[] vhit;
334  // delete[] vt;
335  // return -1;
336  int* LenCltmp = new int[MaxLen];
337  CopyVector(LenCl, LenCltmp, MaxLen);
338  delete[] LenCl;
339  LenCl = new int[MaxLen * 2];
340  ZeroVector(LenCl, MaxLen * 2);
341  CopyVector(LenCltmp, LenCl, MaxLen);
342  delete[] LenCltmp;
343  MaxLen *= 2;
344  // std::cout << "Extend array size to " << MaxLen << std::endl;
345  }
346  nCl++;
347  LenCl[nCl - 1] = next - ib;
348  if (nCl > 1)
349  {
350  // Job to glue the subclusters with common edge
351  //
352  iab = vhit[ib].ich; // The last subcl begin
353  iae = vhit[ie].ich; // The last subcl end
354  last = ib - 1; // The prelast subcl end
355  LastCl = nCl - 2;
356  for (int iCl = LastCl; iCl >= 0; iCl--)
357  {
358  leng = LenCl[iCl];
359 
360  if (iab - vhit[last].ich > fNx) goto new_ich;
361  for (int ichc = last; ichc >= last - leng + 1; ichc--)
362  {
363  // if( iab-vhit[ichc].ich > fNx ) goto new_icl; // From PHENIX version !!! This may be not right for complicated clusters, where tower ordering is not conserved
364 
365  // if( iae-vhit[ichc].ich >= fNx // From PHENIX version
366  if ((vhit[ichc].ich + fNx <= iae && vhit[ichc].ich + fNx >= iab) || (bCYL && (iae % fNx == fNx - 1) && (iae - vhit[ichc].ich == fNx - 1)) // Only for CYLinder geom !!!!
367  )
368  {
369  // Swap iCl-cluster towers (of length "leng") with whatever was between it and the last subcluster (of length "ib-1-last") - to make it adjecent to the last subcluster
370  CopyVector(&vhit[last + 1 - leng], vt, leng);
371  CopyVector(&vhit[last + 1], &vhit[last + 1 - leng], ib - 1 - last);
372  CopyVector(vt, &vhit[ib - leng], leng);
373 
374  // Now the number of clusters is reduced by 1 and the length of the last one increased by iCl-cluster length "leng"
375  for (int i = iCl; i < nCl - 2; i++) LenCl[i] = LenCl[i + 1];
376  ib -= leng;
377  LenCl[nCl - 2] = LenCl[nCl - 1] + leng;
378  nCl--;
379  goto new_icl;
380  }
381  } // for( int ichc=last
382 
383  new_icl:
384  last = last - leng;
385  } // for( int iCl=LastCl
386 
387  } // if( nCl > 1
388 
389  } // if( (ia-vhit
390 
391  new_ich:
392  continue;
393  } // for( ich=1
394 
395  if (nCl > 0)
396  {
397  ib = 0;
398  for (int iCl = 0; iCl < nCl; iCl++)
399  {
400  leng = LenCl[iCl];
401  hl.clear();
402  for (ich = 0; ich < leng; ich++) hl.push_back(vhit[ib + ich]);
403  Clt.ReInitialize(hl);
404  ib += LenCl[iCl];
405  fClusters->push_back(Clt);
406  }
407  }
408  delete[] LenCl;
409  delete[] vhit;
410  delete[] vt;
411 
412  return nCl;
413 }
414 
415 // ///////////////////////////////////////////////////////////////////////////
416 
417 void BEmcRec::Momenta(std::vector<EmcModule>* phit, float& pe, float& px,
418  float& py, float& pxx, float& pyy, float& pyx,
419  float thresh)
420 {
421  // First and second momenta calculation
422 
423  float a, x, y, e, xx, yy, yx;
424  std::vector<EmcModule>::iterator ph;
425 
426  pe = 0;
427  px = 0;
428  py = 0;
429  pxx = 0;
430  pyy = 0;
431  pyx = 0;
432  if (phit->empty()) return;
433 
434  // Find max energy tower
435  //
436  ph = phit->begin();
437  float emax = 0;
438  int ichmax = 0;
439  while (ph != phit->end())
440  {
441  a = ph->amp;
442  if (a > emax)
443  {
444  emax = a;
445  ichmax = ph->ich;
446  }
447  ++ph;
448  }
449  if (emax <= 0) return;
450 
451  int iymax = ichmax / fNx;
452  int ixmax = ichmax - iymax * fNx;
453 
454  // Calculate CG relative to max energy tower
455 
456  x = 0;
457  y = 0;
458  e = 0;
459  xx = 0;
460  yy = 0;
461  yx = 0;
462  ph = phit->begin();
463  while (ph != phit->end())
464  {
465  a = ph->amp;
466  if (a > thresh)
467  {
468  int iy = ph->ich / fNx;
469  int ix = ph->ich - iy * fNx;
470  int idx = iTowerDist(ixmax, ix);
471  int idy = iy - iymax;
472  e += a;
473  x += idx * a;
474  y += idy * a;
475  xx += a * idx * idx;
476  yy += a * idy * idy;
477  yx += a * idx * idy;
478  }
479  ++ph;
480  }
481  pe = e;
482 
483  if (e > 0)
484  {
485  x /= e;
486  y /= e;
487  xx = xx / e - x * x;
488  yy = yy / e - y * y;
489  yx = yx / e - y * x;
490 
491  x += ixmax;
492  y += iymax;
493 
494  while (x < -0.5) x += float(fNx);
495  while (x >= fNx - 0.5) x -= float(fNx);
496 
497  px = x;
498  py = y;
499  pxx = xx;
500  pyy = yy;
501  pyx = yx;
502  }
503 }
504 
505 // ///////////////////////////////////////////////////////////////////////////
506 
507 float BEmcRec::PredictEnergy(float en, float xcg, float ycg, int ix, int iy)
508 {
509  if (_emcprof != nullptr && bProfileProb) return PredictEnergyProb(en, xcg, ycg, ix, iy);
510 
511  float dx = fabs(fTowerDist(float(ix), xcg));
512  float dy = ycg - iy;
513  return PredictEnergyParam(en, dx, dy);
514 }
515 
516 float BEmcRec::PredictEnergyParam(float /*en*/, float xc, float yc)
517 {
518  // Calculates the energy deposited in the tower, the distance between
519  // its center and shower Center of Gravity being (xc,yc)
520  // en - shower energy
521 
522  float dx, dy, r1, r2, r3;
523  float fPpar1, fPpar2, fPpar3, fPpar4;
524 
525  float fPshiftx = 0; // !!!!! Untill tuned ... may not be necessary
526  float fPshifty = 0; // !!!!! Untill tuned ... may not be necessary
527 
528  /*
529  float lgE;
530  if (en <= 1.e-10)
531  lgE = 0;
532  else
533  lgE = log(en);
534  fPpar1=0.59-(1.45+0.13*lgE)*sin2a;
535  fPpar2=0.265+(0.80+0.32*lgE)*sin2a;
536  fPpar3=0.25+(0.45-0.036*lgE)*sin2a;
537  fPpar4=0.42;
538  */
539  fPpar1 = 0.549;
540  fPpar2 = 0.304;
541  fPpar3 = 0.389;
542  fPpar4 = 0.326;
543  /*
544  fPpar1 = 0.486;
545  fPpar2 = 0.302;
546  fPpar3 = 0.354;
547  fPpar4 = 0.407;
548  */
549  /*
550  fPpar1 = 0.343;
551  fPpar2 = 0.509;
552  fPpar3 = 0.199;
553  fPpar4 = 0.548;
554  */
555 
556  // if (en > 0) SetProfileParameters(-1, en, xc, yc);
557 
558  dx = fabs(xc - fPshiftx);
559  dy = fabs(yc - fPshifty);
560  r2 = dx * dx + dy * dy;
561  r1 = sqrt(r2);
562  r3 = r2 * r1;
563  double e = fPpar1 * exp(-r3 / fPpar2) + fPpar3 * exp(-r1 / fPpar4);
564 
565  return e;
566 }
567 
568 float BEmcRec::PredictEnergyProb(float en, float xcg, float ycg, int ix, int iy)
569 // Predict tower energy from profiles used in GetProb()
570 // This is expected to be used in BEmcCluster::GetSubClusters
571 {
572  if (_emcprof == nullptr) return -1;
573 
574  while (xcg < -0.5) xcg += float(fNx);
575  while (xcg >= fNx - 0.5) xcg -= float(fNx);
576 
577  int ixcg = int(xcg + 0.5);
578  int iycg = int(ycg + 0.5);
579  float ddx = fabs(xcg - ixcg);
580  float ddy = fabs(ycg - iycg);
581 
582  float xg, yg, zg;
583  Tower2Global(en, xcg, ycg, xg, yg, zg);
584 
585  float theta, phi;
586  GetImpactThetaPhi(xg, yg, zg, theta, phi);
587 
588  int isx = 1;
589  if (xcg - ixcg < 0) isx = -1;
590  int isy = 1;
591  if (ycg - iycg < 0) isy = -1;
592 
593  int idx = iTowerDist(ixcg, ix) * isx;
594  int idy = (iy - iycg) * isy;
595 
596  int id = -1;
597  if (idx == 0 && idy == 0)
598  id = 0;
599  else if (idx == 1 && idy == 0)
600  id = 1;
601  else if (idx == 1 && idy == 1)
602  id = 2;
603  else if (idx == 0 && idy == 1)
604  id = 3;
605 
606  if (id < 0)
607  {
608  float dx = fabs(fTowerDist(xcg, float(ix)));
609  float dy = fabs(iy - ycg);
610  float rr = sqrt(dx * dx + dy * dy);
611  // return PredictEnergyParam(en, dx, dy);
612  return _emcprof->PredictEnergyR(en, theta, phi, rr);
613  }
614 
615  float ep[4], err[4];
616  for (int ip = 0; ip < 4; ip++)
617  {
618  _emcprof->PredictEnergy(ip, en, theta, phi, ddx, ddy, ep[ip], err[ip]);
619  }
620 
621  float eout;
622 
623  if (id == 0)
624  eout = (ep[1] + ep[2]) / 2. + ep[3];
625  else if (id == 1)
626  eout = (ep[0] - ep[2]) / 2. - ep[3];
627  else if (id == 3)
628  eout = (ep[0] - ep[1]) / 2. - ep[3];
629  else
630  eout = ep[3];
631 
632  // if( eout<0 ) printf("id=%d eout=%f: ep= %f %f %f %f Input: E=%f xcg=%f ycg=%f\n",id,eout,ep[0],ep[1],ep[2],ep[3],en,xcg,ycg);
633  if (eout < 0) eout = 1e-6;
634 
635  return eout;
636 }
637 
638 // ///////////////////////////////////////////////////////////////////////////
639 
640 float BEmcRec::GetTowerEnergy(int iy, int iz, std::vector<EmcModule>* plist)
641 {
642  int nn = plist->size();
643  if (nn <= 0) return 0;
644 
645  for (int i = 0; i < nn; i++)
646  {
647  int ich = (*plist)[i].ich;
648  int iyt = ich / fNx;
649  int izt = ich % fNx;
650  if (iy == iyt && iz == izt)
651  {
652  return (*plist)[i].amp;
653  }
654  }
655  return 0;
656 }
657 
658 // !!!!! Change here to a ponter to HitList
659 float BEmcRec::GetProb(std::vector<EmcModule> HitList, float en, float xg, float yg, float zg, float& chi2, int& ndf)
660 // Do nothing; should be defined in a detector specific module BEmcRec{Name}
661 {
662  // float enoise = 0.01; // 10 MeV per tower
663  float enoise = GetProbNoiseParam();
664  // float thresh = 0.01;
665  float thresh = GetTowerThreshold();
666 
667  chi2 = 0;
668  ndf = 0;
669  if (_emcprof == nullptr) return -1;
670 
671  if (!(_emcprof->IsLoaded()))
672  {
673  return -1;
674  }
675 
676  int nn = HitList.size();
677  if (nn <= 0) return -1;
678 
679  float theta, phi;
680  GetImpactThetaPhi(xg, yg, zg, theta, phi);
681 
682  // z coordinate below means x coordinate
683 
684  float etot;
685  float zcg, ycg;
686  float zz, yy, yz;
687  Momenta(&HitList, etot, zcg, ycg, zz, yy, yz, thresh);
688 
689  int iz0cg = int(zcg + 0.5);
690  int iy0cg = int(ycg + 0.5);
691  float ddz = fabs(zcg - iz0cg);
692  float ddy = fabs(ycg - iy0cg);
693 
694  int isz = 1;
695  if (zcg - iz0cg < 0) isz = -1;
696  int isy = 1;
697  if (ycg - iy0cg < 0) isy = -1;
698 
699  // 4 central towers: 43
700  // 12
701  // Tower 1 - central one
702  float e1, e2, e3, e4;
703  e1 = GetTowerEnergy(iy0cg, iz0cg, &HitList);
704  e2 = GetTowerEnergy(iy0cg, iz0cg + isz, &HitList);
705  e3 = GetTowerEnergy(iy0cg + isy, iz0cg + isz, &HitList);
706  e4 = GetTowerEnergy(iy0cg + isy, iz0cg, &HitList);
707  if (e1 < thresh) e1 = 0;
708  if (e2 < thresh) e2 = 0;
709  if (e3 < thresh) e3 = 0;
710  if (e4 < thresh) e4 = 0;
711 
712  float e1t = (e1 + e2 + e3 + e4) / etot;
713  float e2t = (e1 + e2 - e3 - e4) / etot;
714  float e3t = (e1 - e2 - e3 + e4) / etot;
715  float e4t = (e3) / etot;
716  // float rr = sqrt((0.5-ddz)*(0.5-ddz)+(0.5-ddy)*(0.5-ddy));
717 
718  // Predicted values
719  const int NP = 4; // From BEmcProfile
720  float ep[NP];
721  float err[NP];
722  for (int ip = 0; ip < NP; ip++)
723  {
724  _emcprof->PredictEnergy(ip, en, theta, phi, ddz, ddy, ep[ip], err[ip]);
725  if (ep[ip] < 0)
726  {
727  return -1;
728  }
729  if (ip < 3)
730  {
731  err[ip] = sqrt(err[ip] * err[ip] + 4 * enoise * enoise / etot / etot);
732  }
733  else
734  {
735  err[ip] = sqrt(err[ip] * err[ip] + 1 * enoise * enoise / etot / etot);
736  }
737  }
738 
739  chi2 = 0.;
740  chi2 += (ep[0] - e1t) * (ep[0] - e1t) / err[0] / err[0];
741  chi2 += (ep[1] - e2t) * (ep[1] - e2t) / err[1] / err[1];
742  chi2 += (ep[2] - e3t) * (ep[2] - e3t) / err[2] / err[2];
743  chi2 += (ep[3] - e4t) * (ep[3] - e4t) / err[3] / err[3];
744  ndf = 4;
745 
746  chi2 /= 1.5;
747 
748  float prob = TMath::Prob(chi2, ndf);
749 
750  return prob;
751 }
752 
753 // ///////////////////////////////////////////////////////////////////////////
754 // Static functions
755 
756 int BEmcRec::HitNCompare(const void* h1, const void* h2)
757 {
758  return (static_cast<const EmcModule*>(h1)->ich - static_cast<const EmcModule*>(h2)->ich);
759 }
760 
761 // ///////////////////////////////////////////////////////////////////////////
762 
763 int BEmcRec::HitACompare(const void* h1, const void* h2)
764 {
765  float amp1 = static_cast<const EmcModule*>(h1)->amp;
766  float amp2 = static_cast<const EmcModule*>(h2)->amp;
767  return (amp1 < amp2) ? 1 : (amp1 > amp2) ? -1
768  : 0;
769 }
770 
771 // ///////////////////////////////////////////////////////////////////////////
772 
773 void BEmcRec::ZeroVector(int* v, int N)
774 {
775  int* p = v;
776  for (int i = 0; i < N; i++) *p++ = 0;
777 }
778 
779 // ///////////////////////////////////////////////////////////////////////////
780 
781 void BEmcRec::ZeroVector(float* v, int N)
782 {
783  float* p = v;
784  for (int i = 0; i < N; i++) *p++ = 0;
785 }
786 
787 // ///////////////////////////////////////////////////////////////////////////
788 
790 {
791  // memset(v, 0, N*sizeof(EmcModule));
792  for (int i = 0; i < N; i++)
793  {
794  v[i].ich = 0;
795  v[i].amp = 0;
796  v[i].tof = 0;
797  }
798 }
799 
800 // ///////////////////////////////////////////////////////////////////////////
801 
802 void BEmcRec::CopyVector(int* from, int* to, int N)
803 {
804  if (N <= 0) return;
805  for (int i = 0; i < N; i++) to[i] = from[i];
806 }
807 
808 // ///////////////////////////////////////////////////////////////////////////
809 
811 {
812  if (N <= 0)
813  {
814  return;
815  }
816  for (int i = 0; i < N; i++)
817  {
818  to[i] = from[i];
819  }
820 }
821 
822 // ///////////////////////////////////////////////////////////////////////////
823 
824 /* Future improvements:
825 
826 1. FindClusters(): to ensure that all EmcModules are above energy threshold
827 set by SetThreshold routine (or default one)
828 
829 */
830 
831 // ///////////////////////////////////////////////////////////////////////////
832 // EOF