ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BEmcCluster.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BEmcCluster.cc
1 // Name: EmcCluster.cxx
2 // Author: A. Bazilevsky (RIKEN-BNL)
3 // Major modifications by M. Volkov (RRC KI) Jan 27 2000
4 
5 #include "BEmcCluster.h"
6 #include "BEmcRec.h"
7 
8 #include <iostream>
9 
10 using namespace std;
11 
12 // Define and initialize static members
13 
14 // Max number of peaks in cluster; used in EmcCluster::GetSubClusters(...)
15 int const EmcCluster::fgMaxNofPeaks = 1000;
16 
17 // Used in EmcCluster::GetSubClusters(...), it is the number of iterations
18 // when fitting photons to peak areas
19 int const EmcCluster::fgPeakIter = 6;
20 
21 // Emin cuts cluster energy: Ecl >= Epk1+Epk2 +...+Epkn !!!!!!!!!!!!!
22 float const EmcCluster::fgEmin = 0.002;
23 
25  : ich(0)
26  , amp(0)
27  , tof(0)
28 {
29 }
30 
31 //_____________________________________________________________________________
32 EmcModule::EmcModule(int ich_, float amp_, float tof_)
33  : ich(ich_)
34  , amp(amp_)
35  , tof(tof_)
36 {
37 }
38 
39 // ///////////////////////////////////////////////////////////////////////////
40 // EmcCluster member functions
41 
42 void EmcCluster::GetCorrPos(float& xc, float& yc)
43 // Returns the cluster corrected position in tower units
44 // Corrected for S-oscilations, not for shower depth
45 // Shower depth (z-coord) is defined in fOwner->Tower2Global()
46 {
47  float e, x, y, xx, yy, xy;
48 
49  fOwner->Momenta(&fHitList, e, x, y, xx, yy, xy);
50  fOwner->CorrectPosition(e, x, y, xc, yc);
51 }
52 
53 // ///////////////////////////////////////////////////////////////////////////
54 
55 void EmcCluster::GetGlobalPos(float& xA, float& yA, float& zA)
56 // Returns the cluster position in PHENIX global coord system
57 {
58  float xc, yc;
59  float e = GetTotalEnergy();
60  GetCorrPos(xc, yc);
61  fOwner->Tower2Global(e, xc, yc, xA, yA, zA);
62 }
63 
64 // ///////////////////////////////////////////////////////////////////////////
65 
67 // Returns the energy of the ich-tower (0 if ich not found in the fHitList)
68 {
69  vector<EmcModule>::iterator ph;
70  if (fHitList.empty()) return 0;
71  ph = fHitList.begin();
72  while (ph != fHitList.end())
73  {
74  if ((*ph).ich == ich) return (*ph).amp;
75  ++ph;
76  }
77  return 0;
78 }
79 
80 // ///////////////////////////////////////////////////////////////////////////
81 
82 float EmcCluster::GetTowerEnergy(int ix, int iy)
83 // Returns the energy of the tower ix,iy (0 if tower not found in the fHitList)
84 {
85  vector<EmcModule>::iterator ph;
86 
87  if (fHitList.empty()) return 0;
88  ph = fHitList.begin();
89  int ich = iy * fOwner->GetNx() + ix;
90  return GetTowerEnergy(ich);
91 }
92 
93 // ///////////////////////////////////////////////////////////////////////////
94 
96 // Returns the ToF of the ich-tower (0 if ich not found in the fHitList)
97 {
98  vector<EmcModule>::iterator ph;
99  if (fHitList.empty()) return 0;
100  ph = fHitList.begin();
101  while (ph != fHitList.end())
102  {
103  if ((*ph).ich == ich) return (*ph).tof;
104  ++ph;
105  }
106  return 0;
107 }
108 
109 // ///////////////////////////////////////////////////////////////////////////
110 
112 // Returns the cluster total energy
113 {
114  vector<EmcModule>::iterator ph;
115  float et = 0;
116  if (fHitList.empty()) return 0;
117  ph = fHitList.begin();
118  while (ph != fHitList.end())
119  {
120  et += (*ph).amp;
121  ++ph;
122  }
123  return et;
124 }
125 
126 // ///////////////////////////////////////////////////////////////////////////
127 
129 // Returns the energy in core towers around the cluster Center of Gravity
130 // Corrected for energy leak sidewise from core towers
131 {
132  float e, x, y, xx, yy, xy;
133  float ecore, ecorecorr;
134  ecore = GetECore();
135  fOwner->Momenta(&fHitList, e, x, y, xx, yy, xy);
136  fOwner->CorrectECore(ecore, x, y, ecorecorr);
137  return ecorecorr;
138 }
139 
141 // Returns the energy in core towers around the cluster Center of Gravity
142 {
143  const float thresh = 0.01;
144 
145  vector<EmcModule>::iterator ph;
146  float xcg, ycg, xx, xy, yy;
147  float energy, es;
148 
149  fOwner->Momenta(&fHitList, energy, xcg, ycg, xx, yy, xy);
150  // fOwner->SetProfileParameters(0, energy, xcg, ycg);
151 
152  es = 0;
153  if (fHitList.empty()) return 0;
154  ph = fHitList.begin();
155  while (ph != fHitList.end())
156  {
157  int ixy = (*ph).ich;
158  int iy = ixy / fOwner->GetNx();
159  int ix = ixy - iy * fOwner->GetNx();
160  // dx = xcg - ix;
161  // float dx = fOwner->fTowerDist(float(ix), xcg);
162  // float dy = ycg - iy;
163  // float et = fOwner->PredictEnergy(dx, dy, energy);
164  float et = fOwner->PredictEnergy(energy, xcg, ycg, ix, iy);
165  if (et > thresh) es += (*ph).amp;
166  ++ph;
167  }
168  return es;
169 }
170 
171 // ///////////////////////////////////////////////////////////////////////////
172 
174 // Returns the energy in 2x2 towers around the cluster Center of Gravity
175 {
176  float et, xcg, ycg, xx, xy, yy;
177  float e1, e2, e3, e4;
178  int ix0, iy0, isx, isy;
179 
180  fOwner->Momenta(&fHitList, et, xcg, ycg, xx, yy, xy);
181  ix0 = int(xcg + 0.5);
182  iy0 = int(ycg + 0.5);
183 
184  isx = 1;
185  if (xcg - ix0 < 0) isx = -1;
186  isy = 1;
187  if (ycg - iy0 < 0) isy = -1;
188 
189  e1 = GetTowerEnergy(ix0, iy0);
190  e2 = GetTowerEnergy(ix0 + isx, iy0);
191  e3 = GetTowerEnergy(ix0 + isx, iy0 + isy);
192  e4 = GetTowerEnergy(ix0, iy0 + isy);
193 
194  return (e1 + e2 + e3 + e4);
195 }
196 // ///////////////////////////////////////////////////////////////////////////
197 
199 // Returns the energy in 3x3 towers around the cluster Center of Gravity
200 {
201  float et, xcg, ycg, xx, xy, yy;
202  int ich, ix0, iy0, nhit;
203 
204  nhit = fHitList.size();
205 
206  if (nhit <= 0) return 0;
207 
208  fOwner->Momenta(&fHitList, et, xcg, ycg, xx, yy, xy);
209  ix0 = int(xcg + 0.5);
210  iy0 = int(ycg + 0.5);
211  ich = iy0 * fOwner->GetNx() + ix0;
212 
213  return GetE9(ich);
214 }
215 
216 // ///////////////////////////////////////////////////////////////////////////
217 
218 float EmcCluster::GetE9(int ich)
219 // Returns the energy in 3x3 towers around the tower ich
220 {
221  vector<EmcModule>::iterator ph;
222 
223  int iy0 = ich / fOwner->GetNx();
224  int ix0 = ich - iy0 * fOwner->GetNx();
225 
226  float es = 0;
227  if (fHitList.empty()) return 0;
228  ph = fHitList.begin();
229  while (ph != fHitList.end())
230  {
231  int ixy = (*ph).ich;
232  int iy = ixy / fOwner->GetNx();
233  int ix = ixy - iy * fOwner->GetNx();
234  int idx = fOwner->iTowerDist(ix0, ix);
235  int idy = iy - iy0;
236  if (abs(idx) <= 1 && abs(idy) <= 1) es += (*ph).amp;
237  ++ph;
238  }
239  return es;
240 }
241 
242 // ///////////////////////////////////////////////////////////////////////////
243 
245 // Returns the EmcModule with the maximum energy
246 {
247  vector<EmcModule>::iterator ph;
248  float emax = 0;
249  EmcModule ht;
250 
251  ht.ich = -1;
252  ht.amp = 0;
253  ht.tof = 0;
254  if (fHitList.empty()) return ht;
255 
256  ph = fHitList.begin();
257  while (ph != fHitList.end())
258  {
259  if ((*ph).amp > emax)
260  {
261  emax = (*ph).amp;
262  ht = *ph;
263  }
264  ++ph;
265  }
266  return ht;
267 }
268 
269 // ///////////////////////////////////////////////////////////////////////////
270 
271 void EmcCluster::GetMoments(float& x, float& y, float& xx, float& xy, float& yy)
272 // Returns cluster 1-st (px,py) and 2-d momenta (pxx,pxy,pyy) in cell unit
273 {
274  float e;
275  fOwner->Momenta(&fHitList, e, x, y, xx, yy, xy);
276 }
277 
278 // ///////////////////////////////////////////////////////////////////////////
279 
280 float EmcCluster::GetProb(float& chi2, int& ndf)
281 {
282  float e, xg, yg, zg;
283  e = GetTotalEnergy();
284  GetGlobalPos(xg, yg, zg);
285  return fOwner->GetProb(fHitList, e, xg, yg, zg, chi2, ndf);
286 }
287 
288 // ///////////////////////////////////////////////////////////////////////////
289 
290 int EmcCluster::GetSubClusters(vector<EmcCluster>* PkList, vector<EmcModule>* ppeaks)
291 {
292  // Splits the cluster onto subclusters
293  // The number of subclusters is equal to the number of Local Maxima in a cluster.
294  // Local Maxima can have the energy not less then
295  // defined in fgMinPeakEnergy
296  //
297  // Output: PkList - vector of subclusters
298  // ppeaks - vector of peak EmcModules (one for each subcluster)
299  //
300  // Returns: >= 0 Number of Peaks;
301  // -1 The number of Peaks is greater then fgMaxNofPeaks;
302  // (just increase parameter fgMaxNofPeaks)
303 
304  int npk, ipk, nhit;
305  int ixypk, ixpk, iypk, in, nh, ic;
306  int ixy, ix, iy, nn;
307  int idx, idy;
308  int ig, ng, igmpk1[fgMaxNofPeaks], igmpk2[fgMaxNofPeaks];
309  int PeakCh[fgMaxNofPeaks];
310  float epk[fgMaxNofPeaks * 2], xpk[fgMaxNofPeaks * 2], ypk[fgMaxNofPeaks * 2];
311  float ratio, eg, dx, dy, a;
312  float *Energy[fgMaxNofPeaks], *totEnergy, *tmpEnergy;
313  EmcModule *phit, *hlist, *vv;
314  EmcCluster peak(fOwner);
315  vector<EmcModule>::iterator ph;
316  vector<EmcModule> hl;
317 
318  PkList->clear();
319  ppeaks->clear();
320 
321  nhit = fHitList.size();
322 
323  if (nhit <= 0) return 0;
324 
325  hlist = new EmcModule[nhit];
326 
327  ph = fHitList.begin();
328  vv = hlist;
329  while (ph != fHitList.end()) *vv++ = *ph++;
330 
331  // sort by linear channel number
332  qsort(hlist, nhit, sizeof(EmcModule), fOwner->HitNCompare);
333 
334  //
335  // Find peak (maximum) position (towers with local maximum amp)
336  //
337  npk = 0;
338  for (ic = 0; ic < nhit; ic++)
339  {
340  float amp = hlist[ic].amp;
341  if (amp > fOwner->GetPeakThreshold())
342  {
343  int ich = hlist[ic].ich;
344  //old int ichmax = ich + fOwner->GetNx() + 1;
345  //old int ichmin = ich - fOwner->GetNx() - 1;
346  // Look into three raws only
347  int ichmax = (ich / fOwner->GetNx() + 2) * fOwner->GetNx() - 1;
348  int ichmin = (ich / fOwner->GetNx() - 1) * fOwner->GetNx();
349  int ixc = ich - ich / fOwner->GetNx() * fOwner->GetNx();
350  int in = ic + 1;
351  // check right, and 3 towers above if there is another max
352  while ((in < nhit) && (hlist[in].ich <= ichmax))
353  {
354  int ix = hlist[in].ich - hlist[in].ich / fOwner->GetNx() * fOwner->GetNx();
355  //old if( (abs(ixc-ix) <= 1) && (hlist[in].amp >= amp) ) goto new_ic;
356  if ((abs(fOwner->iTowerDist(ix, ixc)) <= 1) && (hlist[in].amp >= amp)) goto new_ic;
357  in++;
358  }
359  in = ic - 1;
360  // check left, and 3 towers below if there is another max
361  while ((in >= 0) && (hlist[in].ich >= ichmin))
362  {
363  int ix = hlist[in].ich - hlist[in].ich / fOwner->GetNx() * fOwner->GetNx();
364  //old if( (abs(ixc-ix) <= 1) && (hlist[in].amp > amp) ) goto new_ic;
365  if ((abs(fOwner->iTowerDist(ix, ixc)) <= 1) && (hlist[in].amp > amp)) goto new_ic;
366  in--;
367  }
368  if (npk >= fgMaxNofPeaks)
369  {
370  delete[] hlist;
371  cout << "!!! Error in EmcCluster::GetSubClusters(): too many peaks in a cluster (>"
372  << fgMaxNofPeaks
373  << "). May need tower energy threshold increase for clustering." << endl;
374  return -1;
375  }
376 
377  // ic is a maximum in a 3x3 tower group
378  PeakCh[npk] = ic;
379  npk++;
380  }
381  new_ic:
382  continue;
383  }
384  /*
385  for( ipk=0; ipk<npk; ipk++ ) {
386  ic = PeakCh[ipk];
387  cout << " " << ipk << ": E = " << hlist[ic].amp << endl;
388  }
389  */
390 
391  // there was only one peak
392  if (npk <= 1)
393  {
394  hl.clear();
395  for (int ich = 0; ich < nhit; ich++) hl.push_back(hlist[ich]);
396  peak.ReInitialize(hl);
397  PkList->push_back(peak);
398 
399  if (npk == 1)
400  ppeaks->push_back(hlist[PeakCh[0]]);
401  else
402  ppeaks->push_back(GetMaxTower());
403 
404  delete[] hlist;
405  return 1;
406  }
407 
408  // there were more than one peak
409 
410  for (ipk = 0; ipk < npk; ipk++)
411  {
412  Energy[ipk] = new float[nhit];
413  }
414  totEnergy = new float[nhit];
415  tmpEnergy = new float[nhit];
416  for (int i = 0; i < nhit; ++i)
417  {
418  totEnergy[i] = 0.0;
419  tmpEnergy[i] = 0.0;
420  }
421  //
422  // Divide energy in towers among photons positioned in every peak
423  //
424  ratio = 1.;
425  for (int iter = 0; iter < fgPeakIter; iter++)
426  {
427  fOwner->ZeroVector(tmpEnergy, nhit);
428  for (ipk = 0; ipk < npk; ipk++)
429  {
430  ic = PeakCh[ipk];
431  if (iter > 0) ratio = Energy[ipk][ic] / totEnergy[ic];
432  eg = hlist[ic].amp * ratio;
433  ixypk = hlist[ic].ich;
434  iypk = ixypk / fOwner->GetNx();
435  ixpk = ixypk - iypk * fOwner->GetNx();
436  epk[ipk] = eg;
437  xpk[ipk] = 0.; // CoG from the peak tower
438  ypk[ipk] = 0.; // CoG from the peak tower
439 
440  int ichmax = (iypk + 2) * fOwner->GetNx() - 1;
441  int ichmin = (iypk - 1) * fOwner->GetNx();
442 
443  // add up energies of tower to the right and up
444  if (ic < nhit - 1)
445  {
446  for (in = ic + 1; in < nhit; in++)
447  {
448  ixy = hlist[in].ich;
449  iy = ixy / fOwner->GetNx();
450  ix = ixy - iy * fOwner->GetNx();
451 
452  //old if( (ixy-ixypk) > fOwner->GetNx()+1 ) break;
453  if (ixy > ichmax) break;
454  //old if( abs(ix-ixpk) <= 1 ) {
455  idx = fOwner->iTowerDist(ixpk, ix);
456  idy = iy - iypk;
457  if (abs(idx) <= 1)
458  {
459  if (iter > 0) ratio = Energy[ipk][in] / totEnergy[in];
460  eg = hlist[in].amp * ratio;
461  epk[ipk] += eg;
462  xpk[ipk] += eg * idx;
463  ypk[ipk] += eg * idy;
464  }
465  }
466  } // if ic
467 
468  // add up energies of tower to the left and down
469  if (ic >= 1)
470  {
471  for (in = ic - 1; in >= 0; in--)
472  {
473  ixy = hlist[in].ich;
474  iy = ixy / fOwner->GetNx();
475  ix = ixy - iy * fOwner->GetNx();
476 
477  //old if( (ixypk-ixy) > fOwner->GetNx()+1 ) break;
478  if (ixy < ichmin) break;
479  //old if( abs(ix-ixpk) <= 1 ) {
480  idx = fOwner->iTowerDist(ixpk, ix);
481  idy = iy - iypk;
482  if (abs(idx) <= 1)
483  {
484  if (iter > 0) ratio = Energy[ipk][in] / totEnergy[in];
485  eg = hlist[in].amp * ratio;
486  epk[ipk] += eg;
487  xpk[ipk] += eg * idx;
488  ypk[ipk] += eg * idy;
489  }
490  }
491  } // if ic
492 
493  xpk[ipk] = xpk[ipk] / epk[ipk] + ixpk;
494  ypk[ipk] = ypk[ipk] / epk[ipk] + iypk;
495  // fOwner->SetProfileParameters(0, epk[ipk], xpk[ipk], ypk[ipk]);
496 
497  for (in = 0; in < nhit; in++)
498  {
499  ixy = hlist[in].ich;
500  iy = ixy / fOwner->GetNx();
501  ix = ixy - iy * fOwner->GetNx();
502  dx = fOwner->fTowerDist(float(ix), xpk[ipk]);
503  dy = ypk[ipk] - iy;
504  a = 0;
505 
506  // predict energy within 2.5 cell square around local peak
507  if (ABS(dx) < 2.5 && ABS(dy) < 2.5)
508  // a = epk[ipk] * fOwner->PredictEnergy(dx, dy, epk[ipk]);
509  a = epk[ipk] * fOwner->PredictEnergy(epk[ipk], xpk[ipk], ypk[ipk], ix, iy);
510 
511  Energy[ipk][in] = a;
512  tmpEnergy[in] += a;
513  }
514 
515  } // for ipk
516  float flim = 0.000001;
517  for (in = 0; in < nhit; in++)
518  {
519  if (tmpEnergy[in] > flim)
520  {
521  totEnergy[in] = tmpEnergy[in];
522  }
523  else
524  {
525  totEnergy[in] = flim;
526  }
527  }
528  } // for iter
529 
530  phit = new EmcModule[nhit];
531 
532  ng = 0;
533  for (ipk = 0; ipk < npk; ipk++)
534  {
535  nh = 0;
536 
537  // fill phit, the tower distribution for a peak
538  for (in = 0; in < nhit; in++)
539  {
540  if (tmpEnergy[in] > 0)
541  {
542  ixy = hlist[in].ich;
543  a = hlist[in].amp * Energy[ipk][in] / tmpEnergy[in];
544  if (a > fgEmin)
545  {
546  phit[nh].ich = ixy;
547  phit[nh].amp = a;
548  phit[nh].tof = hlist[in].tof; // Not necessary here
549 
550  nh++;
551  }
552  }
553  }
554 
555  // if number hits > 0
556  if (nh > 0)
557  {
558  // !!!!! Exclude Gamma() for now until it is proved working and useful
559  /*
560  chi=fOwner->Chi2Limit(nh)*10;
561  int ndf; // just a plug for changed Gamma parameter list MV 28.01.00
562 
563  fOwner->Gamma(nh, phit,&chi, &chi0, &epk[ng], &xpk[ng], &ypk[ng],
564  &epk[ng+1], &xpk[ng+1], &ypk[ng+1], ndf);
565 
566  igmpk1[ipk]=ng;
567  igmpk2[ipk]=ng;
568  if( epk[ng+1]>0 ) { ng++; igmpk2[ipk]=ng;}
569  */
570  igmpk1[ipk] = ng;
571  igmpk2[ipk] = ng;
572 
573  ng++;
574  }
575  else
576  {
577  igmpk1[ipk] = -1;
578  igmpk2[ipk] = -1;
579  }
580  } // for( ipk=
581 
582  fOwner->ZeroVector(tmpEnergy, nhit);
583  for (ipk = 0; ipk < npk; ipk++)
584  {
585  ig = igmpk1[ipk];
586  if (ig >= 0)
587  {
588  // cout << " " << ipk << ": X = " << xpk[ig]
589  // << " Y = " << ypk[ig] << endl;
590  // fOwner->SetProfileParameters(0, epk[ig], xpk[ig], ypk[ig]);
591  for (in = 0; in < nhit; in++)
592  {
593  Energy[ipk][in] = 0;
594  for (ig = igmpk1[ipk]; ig <= igmpk2[ipk]; ig++)
595  {
596  ixy = hlist[in].ich;
597  iy = ixy / fOwner->GetNx();
598  ix = ixy - iy * fOwner->GetNx();
599  dx = fOwner->fTowerDist(float(ix), xpk[ig]);
600  dy = ypk[ig] - iy;
601  // a = epk[ig] * fOwner->PredictEnergy(dx, dy, epk[ig]);
602  a = epk[ig] * fOwner->PredictEnergy(epk[ig], xpk[ig], ypk[ig], ix, iy);
603  Energy[ipk][in] += a;
604  tmpEnergy[in] += a;
605  }
606  } // for( in
607  } // if( ig >= 0
608  } // for( ipk
609 
610  nn = 0;
611  for (ipk = 0; ipk < npk; ipk++)
612  {
613  nh = 0;
614  for (in = 0; in < nhit; in++)
615  {
616  if (tmpEnergy[in] > 0)
617  {
618  ixy = hlist[in].ich;
619  a = hlist[in].amp * Energy[ipk][in] / tmpEnergy[in];
620  if (a > fgEmin)
621  {
622  phit[nh].ich = ixy;
623  phit[nh].amp = a;
624  phit[nh].tof = hlist[in].tof;
625 
626  nh++;
627  }
628  }
629  } // for( in
630  if (nh > 0)
631  {
632  // *ip++ = hlist[PeakCh[ipk]];
633  ppeaks->push_back(hlist[PeakCh[ipk]]);
634  hl.clear();
635  for (in = 0; in < nh; in++) hl.push_back(phit[in]);
636  peak.ReInitialize(hl);
637  PkList->push_back(peak);
638  nn++;
639  }
640  } // for( ipk
641 
642  delete[] phit;
643  delete[] hlist;
644  for (ipk = 0; ipk < npk; ipk++) delete[] Energy[ipk];
645  delete[] totEnergy;
646  delete[] tmpEnergy;
647 
648  return nn;
649 }