sPhenix code displayed by LXR master/​coresoftware/​offline/​packages/​CaloReco/​BEmcCluster.cc	Source navigation Diff markup Identifier search General search
	Version: master Revert   Change

File indexing completed on 2025-08-06 08:17:31

 // Name: EmcCluster.cxx
 // Author: A. Bazilevsky (RIKEN-BNL)
 // Major modifications by M. Volkov (RRC KI) Jan 27 2000
 
 #include "BEmcCluster.h"
 #include "BEmcRec.h"
 
 #include <iostream>
 
 // Define and initialize static members
 
 // Max number of peaks in cluster; used in EmcCluster::GetSubClusters(...)
 int const EmcCluster::fgMaxNofPeaks = 1000;
 
 // Used in EmcCluster::GetSubClusters(...), it is the number of iterations
 // when fitting photons to peak areas
 int const EmcCluster::fgPeakIter = 6;
 
 // Emin cuts cluster energy: Ecl >= Epk1+Epk2 +...+Epkn !!!!!!!!!!!!!
 float const EmcCluster::fgEmin = 0.002;
 
 // ///////////////////////////////////////////////////////////////////////////
 // EmcCluster member functions
 
 void EmcCluster::GetCorrPos(float& xc, float& yc)
 // Returns the cluster corrected position in tower units
 // Corrected for S-oscilations, not for shower depth
 // Shower depth (z-coord) is defined in fOwner->Tower2Global()
 {
   float e;
   float x;
   float y;
   float xx;
   float yy;
   float xy;
 
   fOwner->Momenta(&fHitList, e, x, y, xx, yy, xy);
   fOwner->CorrectPosition(e, x, y, xc, yc);
 }
 
 // ///////////////////////////////////////////////////////////////////////////
 
 void EmcCluster::GetGlobalPos(float& xA, float& yA, float& zA)
 // Returns the cluster position in PHENIX global coord system
 {
   float xc;
   float yc;
   float e = GetTotalEnergy();
   GetCorrPos(xc, yc);
   fOwner->Tower2Global(e, xc, yc, xA, yA, zA);
 }
 
 // ///////////////////////////////////////////////////////////////////////////
 
 float EmcCluster::GetTowerEnergy(int ich)
 // Returns the energy of the ich-tower (0 if ich not found in the fHitList)
 {
   std::vector<EmcModule>::iterator ph;
   if (fHitList.empty())
   {
     return 0;
   }
   ph = fHitList.begin();
   while (ph != fHitList.end())
   {
     if ((*ph).ich == ich)
     {
       return (*ph).amp;
     }
     ++ph;
   }
   return 0;
 }
 
 // ///////////////////////////////////////////////////////////////////////////
 
 float EmcCluster::GetTowerEnergy(int ix, int iy)
 // Returns the energy of the tower ix,iy (0 if tower not found in the fHitList)
 {
   std::vector<EmcModule>::iterator ph;
 
   if (fHitList.empty())
   {
     return 0;
   }
   ph = fHitList.begin();
   int ich = (iy * fOwner->GetNx()) + ix;
   return GetTowerEnergy(ich);
 }
 
 // ///////////////////////////////////////////////////////////////////////////
 
 float EmcCluster::GetTowerToF(int ich)
 // Returns the ToF of the ich-tower (0 if ich not found in the fHitList)
 {
   std::vector<EmcModule>::iterator ph;
   if (fHitList.empty())
   {
     return 0;
   }
   ph = fHitList.begin();
   while (ph != fHitList.end())
   {
     if ((*ph).ich == ich)
     {
       return (*ph).tof;
     }
     ++ph;
   }
   return 0;
 }
 
 // ///////////////////////////////////////////////////////////////////////////
 
 float EmcCluster::GetTotalEnergy()
 // Returns the cluster total energy
 {
   std::vector<EmcModule>::iterator ph;
   float et = 0;
   if (fHitList.empty())
   {
     return 0;
   }
   ph = fHitList.begin();
   while (ph != fHitList.end())
   {
     et += (*ph).amp;
     ++ph;
   }
   return et;
 }
 
 // ///////////////////////////////////////////////////////////////////////////
 
 float EmcCluster::GetECoreCorrected()
 // Returns the energy in core towers around the cluster Center of Gravity
 // Corrected for energy leak sidewise from core towers
 {
   float e;
   float x;
   float y;
   float xx;
   float yy;
   float xy;
   float ecore;
   float ecorecorr;
   ecore = GetECore();
   fOwner->Momenta(&fHitList, e, x, y, xx, yy, xy);
   fOwner->CorrectECore(ecore, x, y, ecorecorr);
   return ecorecorr;
 }
 
 float EmcCluster::GetECore()
 // Returns the energy in core towers around the cluster Center of Gravity
 {
   const float thresh = 0.01;
 
   std::vector<EmcModule>::iterator ph;
   float xcg;
   float ycg;
   float xx;
   float xy;
   float yy;
   float energy;
   float es;
 
   fOwner->Momenta(&fHitList, energy, xcg, ycg, xx, yy, xy);
   //  fOwner->SetProfileParameters(0, energy, xcg, ycg);
 
   es = 0;
   if (fHitList.empty())
   {
     return 0;
   }
   ph = fHitList.begin();
   while (ph != fHitList.end())
   {
     int ixy = (*ph).ich;
     int iy = ixy / fOwner->GetNx();
     int ix = ixy - (iy * fOwner->GetNx());
     //      dx = xcg - ix;
     //    float dx = fOwner->fTowerDist(float(ix), xcg);
     //    float dy = ycg - iy;
     //    float et = fOwner->PredictEnergy(dx, dy, energy);
     float et = fOwner->PredictEnergy(energy, xcg, ycg, ix, iy);
     if (et > thresh)
     {
       es += (*ph).amp;
     }
     ++ph;
   }
   return es;
 }
 
 // ///////////////////////////////////////////////////////////////////////////
 
 float EmcCluster::GetE4()
 // Returns the energy in 2x2 towers around the cluster Center of Gravity
 {
   float et;
   float xcg;
   float ycg;
   float xx;
   float xy;
   float yy;
   float e1;
   float e2;
   float e3;
   float e4;
   int ix0;
   int iy0;
   int isx;
   int isy;
 
   fOwner->Momenta(&fHitList, et, xcg, ycg, xx, yy, xy);
 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
   ix0 = int(xcg + 0.5);
 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
   iy0 = int(ycg + 0.5);
 
   isx = 1;
   if (xcg - ix0 < 0)
   {
     isx = -1;
   }
   isy = 1;
   if (ycg - iy0 < 0)
   {
     isy = -1;
   }
 
   e1 = GetTowerEnergy(ix0, iy0);
   e2 = GetTowerEnergy(ix0 + isx, iy0);
   e3 = GetTowerEnergy(ix0 + isx, iy0 + isy);
   e4 = GetTowerEnergy(ix0, iy0 + isy);
 
   return (e1 + e2 + e3 + e4);
 }
 // ///////////////////////////////////////////////////////////////////////////
 
 float EmcCluster::GetE9()
 // Returns the energy in 3x3 towers around the cluster Center of Gravity
 {
   float et;
   float xcg;
   float ycg;
   float xx;
   float xy;
   float yy;
   int ich;
   int ix0;
   int iy0;
   int nhit;
 
   nhit = fHitList.size();
 
   if (nhit <= 0)
   {
     return 0;
   }
 
   fOwner->Momenta(&fHitList, et, xcg, ycg, xx, yy, xy);
 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
   ix0 = int(xcg + 0.5);
 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
   iy0 = int(ycg + 0.5);
   ich = iy0 * fOwner->GetNx() + ix0;
 
   return GetE9(ich);
 }
 
 // ///////////////////////////////////////////////////////////////////////////
 
 float EmcCluster::GetE9(int ich)
 // Returns the energy in 3x3 towers around the tower ich
 {
   std::vector<EmcModule>::iterator ph;
 
   int iy0 = ich / fOwner->GetNx();
   int ix0 = ich - (iy0 * fOwner->GetNx());
 
   float es = 0;
   if (fHitList.empty())
   {
     return 0;
   }
   ph = fHitList.begin();
   while (ph != fHitList.end())
   {
     int ixy = (*ph).ich;
     int iy = ixy / fOwner->GetNx();
     int ix = ixy - (iy * fOwner->GetNx());
     int idx = fOwner->iTowerDist(ix0, ix);
     int idy = iy - iy0;
     if (abs(idx) <= 1 && abs(idy) <= 1)
     {
       es += (*ph).amp;
     }
     ++ph;
   }
   return es;
 }
 
 // ///////////////////////////////////////////////////////////////////////////
 
 EmcModule EmcCluster::GetMaxTower()
 // Returns the EmcModule with the maximum energy
 {
   std::vector<EmcModule>::iterator ph;
   float emax = 0;
   EmcModule ht;
 
   ht.ich = -1;
   ht.amp = 0;
   ht.tof = 0;
   if (fHitList.empty())
   {
     return ht;
   }
 
   ph = fHitList.begin();
   while (ph != fHitList.end())
   {
     if ((*ph).amp > emax)
     {
       emax = (*ph).amp;
       ht = *ph;
     }
     ++ph;
   }
   return ht;
 }
 
 // ///////////////////////////////////////////////////////////////////////////
 
 void EmcCluster::GetMoments(float& x, float& y, float& xx, float& xy, float& yy)
 //  Returns cluster 1-st (px,py) and 2-d momenta (pxx,pxy,pyy) in cell unit
 {
   float e;
   fOwner->Momenta(&fHitList, e, x, y, xx, yy, xy);
 }
 
 // ///////////////////////////////////////////////////////////////////////////
 
 float EmcCluster::GetProb(float& chi2, int& ndf)
 {
   float e;
   float xg;
   float yg;
   float zg;
   e = GetTotalEnergy();
   xg = 0; 
   GetGlobalPos(xg, yg, zg);
   return fOwner->GetProb(fHitList, e, xg, yg, zg, chi2, ndf);
 }
 
 // ///////////////////////////////////////////////////////////////////////////
 
 int EmcCluster::GetSubClusters(std::vector<EmcCluster>& PkList, std::vector<EmcModule>& ppeaks, bool dosubclustersplitting)
 {
   // Splits the cluster onto subclusters
   // The number of subclusters is equal to the number of Local Maxima in a cluster.
   // Local Maxima can have the energy not less then
   // defined in fgMinPeakEnergy
   //
   // Output: PkList - vector of subclusters
   //         ppeaks - vector of peak EmcModules (one for each subcluster)
   //
   // Returns: >= 0 Number of Peaks;
   //          -1 The number of Peaks is greater then fgMaxNofPeaks;
   //         (just increase parameter fgMaxNofPeaks)
 
   int npk;
   int ipk;
   int nhit;
   int ixypk;
   int ixpk;
   int iypk;
   int in;
   int nh;
   int ic;
   int ixy;
   int ix;
   int iy;
   int nn;
   int idx;
   int idy;
   int ig;
   int ng;
   int igmpk1[fgMaxNofPeaks];
   int igmpk2[fgMaxNofPeaks];
   int PeakCh[fgMaxNofPeaks];
   float epk[fgMaxNofPeaks * 2];
   float xpk[fgMaxNofPeaks * 2];
   float ypk[fgMaxNofPeaks * 2];
   float ratio;
   float eg;
   float dx;
   float dy;
   float a;
   float *Energy[fgMaxNofPeaks];
   float *totEnergy;
   float *tmpEnergy;
   EmcModule *phit;
   EmcModule *hlist;
   EmcModule *vv;
   EmcCluster peak(fOwner);
   std::vector<EmcModule>::iterator ph;
   std::vector<EmcModule> hl;
 
   PkList.clear();
   ppeaks.clear();
 
   nhit = fHitList.size();
 
   if (nhit <= 0)
   {
     return 0;
   }
 
   hlist = new EmcModule[nhit];
 
   ph = fHitList.begin();
   vv = hlist;
   while (ph != fHitList.end())
   {
     *vv++ = *ph++;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
   }
 
   // sort by linear channel number
   qsort(hlist, nhit, sizeof(EmcModule), BEmcRec::HitNCompare);
 
   //
   //  Find peak (maximum) position (towers with local maximum amp)
   //
 
   int maxc=0;
   float maxamp = 0;
 
   npk = 0;
   for (ic = 0; ic < nhit; ic++)
   {
     float amp = hlist[ic].amp;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
     if (amp > fOwner->GetPeakThreshold())
     {
       int ich = hlist[ic].ich;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
       // old      int ichmax = ich + fOwner->GetNx() + 1;
       // old      int ichmin = ich - fOwner->GetNx() - 1;
       //  Look into three raws only
       int ichmax = ((ich / fOwner->GetNx() + 2) * fOwner->GetNx()) - 1;
       int ichmin = (ich / fOwner->GetNx() - 1) * fOwner->GetNx();
       int ixc = ich - (ich / fOwner->GetNx() * fOwner->GetNx());
       int inA = ic + 1;
       // check right, and 3 towers above if there is another max
       while ((inA < nhit) && (hlist[inA].ich <= ichmax))// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
       {
         int ixA = hlist[inA].ich - (hlist[inA].ich / fOwner->GetNx() * fOwner->GetNx());// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
         // old  if( (abs(ixc-ixA) <= 1) && (hlist[inA].amp >= amp) ) goto new_ic;
         if ((std::abs(fOwner->iTowerDist(ixA, ixc)) <= 1) && (hlist[inA].amp >= amp))// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
         {
           goto new_ic;
         }
         inA++;
       }
       inA = ic - 1;
       // check left, and 3 towers below if there is another max
       while ((inA >= 0) && (hlist[inA].ich >= ichmin))// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
       {
         int ixA = hlist[inA].ich - (hlist[inA].ich / fOwner->GetNx() * fOwner->GetNx());// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
         // old  if( (abs(ixc-ixA) <= 1) && (hlist[inA].amp > amp) ) goto new_ic;
         if ((std::abs(fOwner->iTowerDist(ixA, ixc)) <= 1) && (hlist[inA].amp > amp))// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
         {
           goto new_ic;
         }
         inA--;
       }
       if (npk >= fgMaxNofPeaks)
       {
         delete[] hlist;
         std::cout << "!!! Error in EmcCluster::GetSubClusters(): too many peaks in a cluster (>"
                   << fgMaxNofPeaks
                   << "). May need tower energy threshold increase for clustering." << std::endl;
         return -1;
       }
 
       // ic is a maximum in a 3x3 tower group
       if (amp > maxamp) {
         maxamp = amp;
         maxc = npk;
       }   
       PeakCh[npk] = ic;
       npk++;
     }
   new_ic:
     continue;
   }
   /*
   for( ipk=0; ipk<npk; ipk++ ) {
     ic = PeakCh[ipk];
     std::cout << "  " << ipk << ": E = " << hlist[ic].amp << std::endl;
   }
   */
 
   if(!dosubclustersplitting){
     hl.clear();
     for (int ich = 0; ich < nhit; ich++)
     {
       hl.push_back(hlist[ich]);// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
     }
     peak.ReInitialize(hl);
     PkList.push_back(peak);
 
     if (npk > 0)
     {
       ppeaks.push_back(hlist[PeakCh[maxc]]);// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
     }
     else
     {
       ppeaks.push_back(GetMaxTower());
     }
     delete[] hlist;
     return 1;
   }
 
 
   // there was only one peak
   if(npk <= 1)
   {
     hl.clear();
     for (int ich = 0; ich < nhit; ich++)
     {
       hl.push_back(hlist[ich]);// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
     }
     peak.ReInitialize(hl);
     PkList.push_back(peak);
 
     if (npk == 1)
     {
       ppeaks.push_back(hlist[PeakCh[0]]);// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
     }
     else
     {
       ppeaks.push_back(GetMaxTower());
     }
 
     delete[] hlist;
     return 1;
   }
 
   // there were more than one peak
 
   for (ipk = 0; ipk < npk; ipk++)
   {
     Energy[ipk] = new float[nhit];
   }
   totEnergy = new float[nhit];
   tmpEnergy = new float[nhit];
   for (int i = 0; i < nhit; ++i)
   {
     totEnergy[i] = 0.0;
     tmpEnergy[i] = 0.0;
   }
   //
   // Divide energy in towers among photons positioned in every peak
   //
   ratio = 1.;
   for (int iter = 0; iter < fgPeakIter; iter++)
   {
     BEmcRec::ZeroVector(tmpEnergy, nhit);
     for (ipk = 0; ipk < npk; ipk++)
     {
       ic = PeakCh[ipk];
       if (iter > 0)
       {
         ratio = Energy[ipk][ic] / totEnergy[ic];
       }
       eg = hlist[ic].amp * ratio;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
       ixypk = hlist[ic].ich;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
       iypk = ixypk / fOwner->GetNx();
       ixpk = ixypk - iypk * fOwner->GetNx();
       epk[ipk] = eg;
       xpk[ipk] = 0.;  // CoG from the peak tower
       ypk[ipk] = 0.;  // CoG from the peak tower
 
       int ichmax = ((iypk + 2) * fOwner->GetNx()) - 1;
       int ichmin = (iypk - 1) * fOwner->GetNx();
 
       // add up energies of tower to the right and up
       if (ic < nhit - 1)
       {
         for (in = ic + 1; in < nhit; in++)
         {
           ixy = hlist[in].ich;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
           iy = ixy / fOwner->GetNx();
           ix = ixy - iy * fOwner->GetNx();
 
           // old      if( (ixy-ixypk) > fOwner->GetNx()+1 ) break;
           if (ixy > ichmax)
           {
             break;
           }
           // old      if( abs(ix-ixpk) <= 1 ) {
           idx = fOwner->iTowerDist(ixpk, ix);
           idy = iy - iypk;
           if (abs(idx) <= 1)
           {
             if (iter > 0)
             {
               ratio = Energy[ipk][in] / totEnergy[in];
             }
             eg = hlist[in].amp * ratio;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
             epk[ipk] += eg;
             xpk[ipk] += eg * idx;
             ypk[ipk] += eg * idy;
           }
         }
       }  // if ic
 
       // add up energies of tower to the left and down
       if (ic >= 1)
       {
         for (in = ic - 1; in >= 0; in--)
         {
           ixy = hlist[in].ich;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
           iy = ixy / fOwner->GetNx();
           ix = ixy - iy * fOwner->GetNx();
 
           // old      if( (ixypk-ixy) > fOwner->GetNx()+1 ) break;
           if (ixy < ichmin)
           {
             break;
           }
           // old      if( abs(ix-ixpk) <= 1 ) {
           idx = fOwner->iTowerDist(ixpk, ix);
           idy = iy - iypk;
           if (abs(idx) <= 1)
           {
             if (iter > 0)
             {
               ratio = Energy[ipk][in] / totEnergy[in];
             }
             eg = hlist[in].amp * ratio;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
             epk[ipk] += eg;
             xpk[ipk] += eg * idx;
             ypk[ipk] += eg * idy;
           }
         }
       }  // if ic
 
       xpk[ipk] = xpk[ipk] / epk[ipk] + ixpk;
       ypk[ipk] = ypk[ipk] / epk[ipk] + iypk;
       //      fOwner->SetProfileParameters(0, epk[ipk], xpk[ipk], ypk[ipk]);
 
       for (in = 0; in < nhit; in++)
       {
         ixy = hlist[in].ich;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
         iy = ixy / fOwner->GetNx();
         ix = ixy - iy * fOwner->GetNx();
         dx = fOwner->fTowerDist(float(ix), xpk[ipk]);
         dy = ypk[ipk] - iy;
         a = 0;
 
         // predict energy within 2.5 cell square around local peak
         if (ABS(dx) < 2.5 && ABS(dy) < 2.5)
         {
           //          a = epk[ipk] * fOwner->PredictEnergy(dx, dy, epk[ipk]);
           a = epk[ipk] * fOwner->PredictEnergy(epk[ipk], xpk[ipk], ypk[ipk], ix, iy);
         }
 
         Energy[ipk][in] = a;
         tmpEnergy[in] += a;
       }
 
     }  // for ipk
     float flim = 0.000001;
     for (in = 0; in < nhit; in++)
     {
       if (tmpEnergy[in] > flim)
       {
         totEnergy[in] = tmpEnergy[in];
       }
       else
       {
         totEnergy[in] = flim;
       }
     }
   }  // for iter
 
   phit = new EmcModule[nhit];
 
   ng = 0;
   for (ipk = 0; ipk < npk; ipk++)
   {
     nh = 0;
 
     // fill phit, the tower distribution for a peak
     for (in = 0; in < nhit; in++)
     {
       if (tmpEnergy[in] > 0)
       {
         ixy = hlist[in].ich;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
         a = hlist[in].amp * Energy[ipk][in] / tmpEnergy[in];// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
         if (a > fgEmin)
         {
           phit[nh].ich = ixy;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
           phit[nh].amp = a;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
           phit[nh].tof = hlist[in].tof; // NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object) // Not necessary here
 
           nh++;
         }
       }
     }
 
     // if number hits > 0
     if (nh > 0)
     {
       // !!!!! Exclude Gamma() for now until it is proved working and useful
       /*
       chi=fOwner->Chi2Limit(nh)*10;
       int ndf; // just a plug for changed Gamma parameter list MV 28.01.00
 
       fOwner->Gamma(nh, phit,&chi, &chi0, &epk[ng], &xpk[ng], &ypk[ng],
                     &epk[ng+1], &xpk[ng+1], &ypk[ng+1], ndf);
 
       igmpk1[ipk]=ng;
       igmpk2[ipk]=ng;
       if( epk[ng+1]>0 ) { ng++; igmpk2[ipk]=ng;}
       */
       igmpk1[ipk] = ng;
       igmpk2[ipk] = ng;
 
       ng++;
     }
     else
     {
       igmpk1[ipk] = -1;
       igmpk2[ipk] = -1;
     }
   }  // for( ipk=
 
   BEmcRec::ZeroVector(tmpEnergy, nhit);
   for (ipk = 0; ipk < npk; ipk++)
   {
     ig = igmpk1[ipk];
     if (ig >= 0)
     {
       //      std::cout << "  " << ipk << ": X = " << xpk[ig]
       //           << " Y = " << ypk[ig] << std::endl;
       //      fOwner->SetProfileParameters(0, epk[ig], xpk[ig], ypk[ig]);
       for (in = 0; in < nhit; in++)
       {
         Energy[ipk][in] = 0;
         for (ig = igmpk1[ipk]; ig <= igmpk2[ipk]; ig++)
         {
           ixy = hlist[in].ich;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
           iy = ixy / fOwner->GetNx();
           ix = ixy - iy * fOwner->GetNx();
           dx = fOwner->fTowerDist(float(ix), xpk[ig]);
           dy = ypk[ig] - iy;
           //          a = epk[ig] * fOwner->PredictEnergy(dx, dy, epk[ig]);
           a = epk[ig] * fOwner->PredictEnergy(epk[ig], xpk[ig], ypk[ig], ix, iy);
           Energy[ipk][in] += a;
           tmpEnergy[in] += a;
         }
       }  // for( in
     }    // if( ig >= 0
   }      // for( ipk
 
   nn = 0;
   for (ipk = 0; ipk < npk; ipk++)
   {
     nh = 0;
     for (in = 0; in < nhit; in++)
     {
       if (tmpEnergy[in] > 0)
       {
         ixy = hlist[in].ich;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
         a = hlist[in].amp * Energy[ipk][in] / tmpEnergy[in];// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
         if (a > fgEmin)
         {
           phit[nh].ich = ixy;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
           phit[nh].amp = a;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
           phit[nh].tof = hlist[in].tof;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
 
           nh++;
         }
       }
     }  // for( in
     if (nh > 0)
     {
       //      *ip++ = hlist[PeakCh[ipk]];
       ppeaks.push_back(hlist[PeakCh[ipk]]);// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
       hl.clear();
       for (in = 0; in < nh; in++)
       {
         hl.push_back(phit[in]);// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
       }
       peak.ReInitialize(hl);
       PkList.push_back(peak);
       nn++;
     }
   }  // for( ipk
 
   delete[] phit;
   delete[] hlist;
   for (ipk = 0; ipk < npk; ipk++)
   {
     delete[] Energy[ipk];
   }
   delete[] totEnergy;
   delete[] tmpEnergy;
 
   return nn;
 }

This page was automatically generated by the 2.3.7 LXR engine. The LXR team