packages/CaloReco/BEmcRec.cc

0001 // Name: BEmcRec.h
0002 // Author: A. Bazilevsky, Apr 2012
0003 // Modified from EmcSectorRec.cxx and EmcScSectorRec.cxx
0004
0005 // BEmcRec -- base class for sPHENIX EMCal
0006
0007 // ///////////////////////////////////////////////////////////////////////////
0008
0009 #include "BEmcRec.h"
0010 #include "BEmcCluster.h"
0011 #include "BEmcProfile.h"
0012
0013 #include <TMath.h>
0014
0015 #include <cmath>
0016 #include <cstdlib>
0017 #include <fstream>
0018 #include <iostream>
0019 #include <utility>
0020
0021 // ///////////////////////////////////////////////////////////////////////////
0022 // BEmcRec member functions
0023
0024 BEmcRec::BEmcRec() : fModules(new std::vector<EmcModule>), fClusters(new std::vector<EmcCluster>)
0025 {
0026   fTowerGeom.clear();
0027
0028
0029 }
0030
0031 // ///////////////////////////////////////////////////////////////////////////
0032
0033 BEmcRec::~BEmcRec()
0034 {
0035   if (fModules)
0036   {
0037     fModules->clear();
0038     delete fModules;
0039   }
0040
0041   if (fClusters)
0042   {
0043     fClusters->clear();
0044     delete fClusters;
0045   }
0046
0047   delete _emcprof;
0048 }
0049
0050 // ///////////////////////////////////////////////////////////////////////////
0051
0052 void BEmcRec::LoadProfile(const std::string& /*fname*/)
0053 {
0054   std::cout << "Warning from BEmcRec::LoadProfile(): No acton defined for shower profile evaluation; should be defined in a detector specific module " << Name() << std::endl;
0055 }
0056
0057 void BEmcRec::PrintTowerGeometry(const std::string& fname)
0058 {
0059   std::ofstream outfile(fname);
0060   if (!outfile.is_open())
0061   {
0062     std::cout << "Error in BEmcRec::PrintTowerGeometry(): Failed to open file "
0063               << fname << std::endl;
0064     return;
0065   }
0066   outfile << "Number of bins:" << std::endl;
0067   outfile << fNx << " " << fNy << std::endl;
0068   outfile << "ix iy x y z dx0 dy0 dz0 dx1 dy1 dz1" << std::endl;
0069   int ich;
0070   TowerGeom geom;
0071   std::map<int, TowerGeom>::iterator it;
0072   for (int iy = 0; iy < fNy; iy++)
0073   {
0074     for (int ix = 0; ix < fNx; ix++)
0075     {
0076       ich = iy * fNx + ix;
0077       it = fTowerGeom.find(ich);
0078       if (it != fTowerGeom.end())
0079       {
0080         geom = it->second;
0081         outfile << ix << " " << iy << " " << geom.Xcenter << " "
0082                 << geom.Ycenter << " " << geom.Zcenter << " " << geom.dX[0] << " "
0083                 << geom.dY[0] << " " << geom.dZ[0] << " " << geom.dX[1] << " "
0084                 << geom.dY[1] << " " << geom.dZ[1] << std::endl;
0085         //  std::cout << "Z0: " << geom.dZ[0] << " || Z1: " << geom.dZ[1] << std::endl;
0086       }
0087     }
0088   }
0089 }
0090
0091 void BEmcRec::PrintTowerGeometryDetailed(const std::string& fname)
0092 {
0093   std::ofstream outfile(fname);
0094   if (!outfile.is_open())
0095   {
0096     std::cout << "Error in BEmcRec::PrintTowerGeometry(): Failed to open file "
0097               << fname << std::endl;
0098     return;
0099   }
0100   outfile << "Number of bins:" << std::endl;
0101   outfile << fNx << " " << fNy << std::endl;
0102   outfile << "ieta iphi vtx_1_x vtx_1_y vtx_1_z vtx_2_x vtx_2_y vtx_2_z vtx_3_x vtx_3_y vtx_3_z vtx_4_x vtx_4_y vtx_4_z vtx_5_x vtx_5_y vtx_5_z vtx_6_x vtx_6_y vtx_6_z vtx_7_x vtx_7_y vtx_7_z vtx_8_x vtx_8_y vtx_8_z\n";
0103   int ich;
0104   RawTowerGeom *geom;
0105   std::map<int, RawTowerGeom*>::iterator it;
0106   for (int iy = 0; iy < fNy; iy++)
0107   {
0108     for (int ix = 0; ix < fNx; ix++)
0109     {
0110       ich = iy * fNx + ix;
0111       it = fTowerGeomDetailed.find(ich);
0112       if (it != fTowerGeomDetailed.end())
0113       {
0114         geom = it->second;
0115         outfile << ix << " " << iy << " ";
0116         for (int ivtx = 0; ivtx < 8; ivtx++) {
0117           outfile << geom->get_vertex_x(ivtx) << " " << geom->get_vertex_y(ivtx) << " " << geom->get_vertex_z(ivtx) << " ";
0118         }
0119         outfile << "\n";
0120       }
0121     }
0122   }
0123 }
0124
0125 bool BEmcRec::GetTowerGeometry(int ix, int iy, TowerGeom& geom)
0126 {
0127   if (ix < 0 || ix >= fNx || iy < 0 || iy >= fNy)
0128   {
0129     return false;
0130   }
0131
0132   int ich = (iy * fNx) + ix;
0133   std::map<int, TowerGeom>::iterator it = fTowerGeom.find(ich);
0134   if (it == fTowerGeom.end())
0135   {
0136     return false;
0137   }
0138
0139   geom = it->second;
0140   return true;
0141 }
0142
0143 bool BEmcRec::SetTowerGeometry(int ix, int iy, const RawTowerGeom& raw_geom0)
0144 {
0145   if (ix < 0 || ix >= fNx || iy < 0 || iy >= fNy)
0146   {
0147     return false;
0148   }
0149
0150   TowerGeom geom; // (intermediate geometry used for S-corrections)
0151   geom.Xcenter = raw_geom0.get_center_x();
0152   geom.Ycenter = raw_geom0.get_center_y();
0153   geom.Zcenter = raw_geom0.get_center_z();
0154
0155   if (m_UseDetailedGeometry)
0156   {
0157     geom.rotX = raw_geom0.get_rotx();
0158     geom.rotY = raw_geom0.get_roty();
0159     geom.rotZ = raw_geom0.get_rotz();
0160
0161     // Describe the eta and phi direction within the tower
0162     geom.dX[0] = raw_geom0.get_center_high_phi_x() - raw_geom0.get_center_low_phi_x();
0163     geom.dY[0] = raw_geom0.get_center_high_phi_y() - raw_geom0.get_center_low_phi_y();
0164     geom.dZ[0] = raw_geom0.get_center_high_phi_z() - raw_geom0.get_center_low_phi_z();
0165     geom.dX[1] = raw_geom0.get_center_high_eta_x() - raw_geom0.get_center_low_eta_x();
0166     geom.dY[1] = raw_geom0.get_center_high_eta_y() - raw_geom0.get_center_low_eta_y();
0167     geom.dZ[1] = raw_geom0.get_center_high_eta_z() - raw_geom0.get_center_low_eta_z();
0168   }
0169   else
0170   {
0171     // By default, an approximate value for the tower rotation will be given
0172     geom.rotX = 0;
0173     geom.rotY = 0;
0174     geom.rotZ = 0;
0175
0176     // The tower eta and phi directions should be computed
0177     // in the CompleteTowerGeometry method afterwards
0178     geom.dX[0] = geom.dX[1] = 0;
0179     geom.dY[0] = geom.dY[1] = 0;
0180     geom.dZ[0] = geom.dZ[1] = 0;
0181   }
0182
0183   int ich = (iy * fNx) + ix;
0184   fTowerGeom[ich] = geom;
0185
0186   if (m_UseDetailedGeometry)
0187   {
0188     // copy the input geometry inside fTowerGeomDetailed (only for debugging)
0189     RawTowerGeomv5 *geom_detailed = new RawTowerGeomv5(raw_geom0);
0190     fTowerGeomDetailed[ich] = geom_detailed;
0191   }
0192
0193   return true;
0194 }
0195
0196 bool BEmcRec::CompleteTowerGeometry()
0197 // Calculates tower front size from coordinates of tower center coordinates
0198 {
0199   if (fTowerGeom.empty() || fNx <= 0)
0200   {
0201     std::cout << "Error in BEmcRec::CalculateTowerSize(): Tower geometry not well setup (NX = "
0202               << fNx << ")" << std::endl;
0203     return false;
0204   }
0205
0206   const int nb = 8;
0207   int idx[nb] = {0, 1, 0, -1, -1, 1, 1, -1};
0208   int idy[nb] = {-1, 0, 1, 0, -1, -1, 1, 1};
0209
0210   std::map<int, TowerGeom>::iterator it;
0211
0212   for (it = fTowerGeom.begin(); it != fTowerGeom.end(); ++it)
0213   {
0214     int ich = it->first;
0215     TowerGeom geom0 = it->second;
0216     int ix = ich % fNx;
0217     int iy = ich / fNx;
0218
0219     TowerGeom geomx;
0220     int inx = 0;
0221
0222     while (inx < nb && (idx[inx] == 0 || !GetTowerGeometry(ix + idx[inx], iy + idy[inx], geomx)))
0223     {
0224       inx++;
0225     }
0226     if (inx >= nb)
0227     {
0228       std::cout << "Error in BEmcRec::CompleteTowerGeometry(): Error when locating neighbour for (ix,iy)=("
0229                 << ix << "," << iy << ")" << std::endl;
0230       return false;
0231     }
0232
0233     TowerGeom geomy;
0234     int iny = 0;
0235
0236     while (iny < nb && (idy[iny] == 0 || !GetTowerGeometry(ix + idx[iny], iy + idy[iny], geomy)))
0237     {
0238       iny++;
0239     }
0240     if (iny >= nb)
0241     {
0242       std::cout << "Error in BEmcRec::CompleteTowerGeometry(): Error when locating neighbour for (ix,iy)=("
0243                 << ix << "," << iy << ")" << std::endl;
0244       return false;
0245     }
0246
0247     geom0.dX[0] = (geomx.Xcenter - geom0.Xcenter) / float(idx[inx]);
0248     geom0.dY[0] = (geomx.Ycenter - geom0.Ycenter) / float(idx[inx]);
0249     geom0.dZ[0] = (geomx.Zcenter - geom0.Zcenter) / float(idx[inx]);
0250     geom0.dX[1] = (geomy.Xcenter - geom0.Xcenter) / float(idy[iny]);
0251     geom0.dY[1] = (geomy.Ycenter - geom0.Ycenter) / float(idy[iny]);
0252     geom0.dZ[1] = (geomy.Zcenter - geom0.Zcenter) / float(idy[iny]);
0253
0254     it->second = geom0;
0255
0256   }  // it = fTowerGeom.begin()
0257
0258   return true;
0259 }
0260
0261 void BEmcRec::Tower2Global(float E, float xC, float yC,
0262                            float& xA, float& yA, float& zA)
0263 // xC and yC are local position in tower units
0264 // For CYL geometry (xC,yC) is actually (phiC,zC)
0265 {
0266
0267
0268   bool flagDoSD =  true;
0269   if (xA < -999)
0270   {
0271     flagDoSD = false;
0272   }
0273
0274   xA = 0;
0275   yA = 0;
0276   zA = 0;
0277
0278 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
0279   int ix = xC + 0.5;  // tower #
0280   if (ix < 0 || ix >= fNx)
0281   {
0282     std::cout << m_ThisName << " Error in BEmcRec::Tower2Global: wrong input x: " << ix << std::endl;
0283     return;
0284   }
0285
0286 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
0287   int iy = yC + 0.5;  // tower #
0288   if (iy < 0 || iy >= fNy)
0289   {
0290     std::cout << "Error in BEmcRec::Tower2Global: wrong input y: " << iy << std::endl;
0291     return;
0292   }
0293
0294   // Get tower where the shower is positioned
0295   TowerGeom geom0;
0296
0297   if (!GetTowerGeometry(ix, iy, geom0))
0298   {
0299     // Weird case: cluster center of gravity outside the EMCal, take geometry from the neighbouring tower
0300     const int idx[4] = {1, 0, -1, 0};
0301     const int idy[4] = {0, 1, 0, -1};
0302     int ii = 0;
0303     while (ii < 4 && !GetTowerGeometry(ix + idx[ii], iy + idy[ii], geom0))
0304     {
0305       ii++;
0306     }
0307     if (ii >= 4)
0308     {
0309       std::cout << "Error in BEmcRec::Tower2Global: can not identify neighbour for tower ("
0310                 << ix << "," << iy << ")" << std::endl;
0311       return;
0312     }
0313     float Xc = geom0.Xcenter - (idx[ii] * geom0.dX[0]) - (idy[ii] * geom0.dX[1]);
0314     float Yc = geom0.Ycenter - (idx[ii] * geom0.dY[0]) - (idy[ii] * geom0.dY[1]);
0315     float Zc = geom0.Zcenter - (idx[ii] * geom0.dZ[0]) - (idy[ii] * geom0.dZ[1]);
0316     geom0.Xcenter = Xc;
0317     geom0.Ycenter = Yc;
0318     geom0.Zcenter = Zc;
0319   }
0320
0321   float xt = geom0.Xcenter + ((xC - ix) * geom0.dX[0]) + ((yC - iy) * geom0.dX[1]);
0322   float yt = geom0.Ycenter + ((xC - ix) * geom0.dY[0]) + ((yC - iy) * geom0.dY[1]);
0323   float zt = geom0.Zcenter + ((xC - ix) * geom0.dZ[0]) + ((yC - iy) * geom0.dZ[1]);
0324
0325   if (flagDoSD)
0326   {
0327     CorrectShowerDepth(ix, iy, E, xt, yt, zt, xA, yA, zA);
0328   }
0329   else
0330   {
0331     // If zvtx is not used (m_UseAltZvtx = 2 in RawClusterBuilderTemplate)
0332     // then there is no correction of z_cluster
0333     float savzt = zt;
0334     CorrectShowerDepth(ix, iy, E, xt, yt, zt, xA, yA, zA);
0335     zA = savzt;
0336   }
0337 }
0338
0339 // ///////////////////////////////////////////////////////////////////////////
0340
0341 int BEmcRec::iTowerDist(int ix1, int ix2) const
0342 // Distrance in tower units
0343 {
0344   int idist = ix2 - ix1;
0345   if (bCYL)
0346   {
0347     int idistr = fNx - abs(idist);  // Always >0
0348     if (idistr < abs(idist))
0349     {  // Then count in opposite direction
0350       if (idist < 0)
0351       {
0352         idist = idistr;
0353       }
0354       else
0355       {
0356         idist = -idistr;
0357       }
0358     }
0359   }
0360   //  std::cout << "Dist " << ix1 << " " << ix2 << ": " << idist << std::endl;
0361   return idist;
0362 }
0363
0364 float BEmcRec::fTowerDist(float x1, float x2) const
0365 {
0366   float dist = x2 - x1;
0367   if (bCYL)
0368   {
0369     float distr = fNx - std::fabs(dist);  // Always >0
0370     if (distr < std::abs(dist))
0371     {  // Then count in opposite direction
0372       if (dist < 0)
0373       {
0374         dist = distr;
0375       }
0376       else
0377       {
0378         dist = -distr;
0379       }
0380     }
0381   }
0382   return dist;
0383 }
0384
0385 // ///////////////////////////////////
0386
0387 int BEmcRec::FindClusters()
0388 // Cluster search algorithm based on Lednev's one developed for GAMS.
0389 // Returns number of clusters found
0390 {
0391   int nhit;
0392   int nCl;
0393   //  int LenCl[fgMaxLen];
0394   int* LenCl;
0395   int next;
0396   int ib;
0397   int ie;
0398   int iab;
0399   int iae;
0400   int last;
0401   int LastCl;
0402   int leng;
0403   int ich;
0404   int ia = 0;
0405
0406   EmcModule* vv;
0407   EmcModule *vhit;
0408   EmcModule *vt;
0409   EmcCluster Clt(this);
0410   std::vector<EmcModule>::iterator ph;
0411   std::vector<EmcModule> hl;
0412
0413   (*fClusters).clear();
0414   nhit = (*fModules).size();
0415
0416   if (nhit <= 0)
0417   {
0418     return 0;
0419   }
0420   if (nhit == 1)
0421   {
0422     Clt.ReInitialize((*fModules));
0423     fClusters->push_back(Clt);
0424     return 1;
0425   }
0426
0427   int MaxLen = fgMaxLen;
0428   LenCl = new int[MaxLen];
0429   ZeroVector(LenCl, MaxLen);
0430
0431   vt = new EmcModule[nhit];
0432   vhit = new EmcModule[nhit];
0433
0434   ph = (*fModules).begin();
0435   vv = vhit;
0436   while (ph != (*fModules).end())
0437   {
0438     *vv++ = *ph++;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
0439   }
0440
0441   qsort(vhit, nhit, sizeof(EmcModule), HitNCompare);
0442
0443   nCl = 0;
0444   next = 0;
0445   for (ich = 1; ich < nhit + 1; ich++)
0446   {
0447     if (ich < nhit)
0448     {
0449       ia = vhit[ich].ich; // NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
0450     }
0451
0452     // New subcluster
0453     // NOLINTNEXTLINE(bugprone-pointer-arithmetic-on-polymorphic-object)
0454     if ((ia - vhit[ich - 1].ich > 1)  // the beginning of new subcluster
0455         || (ich >= nhit)              // just finish defining last sub-cluster
0456         || (ia - ia / fNx * fNx == 0))
0457     {  // new raw -> new subcluster
0458
0459       ib = next;
0460       ie = ich - 1;
0461       next = ich;
0462       if (nCl >= MaxLen)
0463       {
0464         //        delete[] vhit;
0465         //        delete[] vt;
0466         //        return -1;
0467         int* LenCltmp = new int[MaxLen];
0468         CopyVector(LenCl, LenCltmp, MaxLen);
0469         delete[] LenCl;
0470 // NOLINTNEXTLINE(bugprone-implicit-widening-of-multiplication-result)
0471         LenCl = new int[MaxLen * 2];
0472         ZeroVector(LenCl, MaxLen * 2);
0473         CopyVector(LenCltmp, LenCl, MaxLen);
0474         delete[] LenCltmp;
0475         MaxLen *= 2;
0476         //  std::cout << "Extend array size to " << MaxLen << std::endl;
0477       }
0478       nCl++;
0479       LenCl[nCl - 1] = next - ib;
0480       if (nCl > 1)
0481       {
0482         // Job to glue the subclusters with common edge
0483         //
0484         iab = vhit[ib].ich;  // NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object) // The last subcl begin
0485         iae = vhit[ie].ich;  // NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object) // The last subcl end
0486         last = ib - 1;       // The prelast subcl end
0487         LastCl = nCl - 2;
0488         for (int iCl = LastCl; iCl >= 0; iCl--)
0489         {
0490           leng = LenCl[iCl];
0491
0492           if (iab - vhit[last].ich > fNx)// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
0493           {
0494             goto new_ich;
0495           }
0496           for (int ichc = last; ichc >= last - leng + 1; ichc--)
0497           {
0498             //      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
0499
0500             //      if( iae-vhit[ichc].ich >= fNx // From PHENIX version
0501         // NOLINTNEXTLINE(bugprone-pointer-arithmetic-on-polymorphic-object)
0502             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 !!!!
0503             )
0504             {
0505               // 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
0506               CopyVector(&vhit[last + 1 - leng], vt, leng);// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
0507               CopyVector(&vhit[last + 1], &vhit[last + 1 - leng], ib - 1 - last);// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
0508               CopyVector(vt, &vhit[ib - leng], leng);// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
0509
0510               // Now the number of clusters is reduced by 1 and the length of the last one increased by iCl-cluster length "leng"
0511               for (int i = iCl; i < nCl - 2; i++)
0512               {
0513                 LenCl[i] = LenCl[i + 1];
0514               }
0515               ib -= leng;
0516               LenCl[nCl - 2] = LenCl[nCl - 1] + leng;
0517               nCl--;
0518               goto new_icl;
0519             }
0520           }  // for( int ichc=last
0521
0522         new_icl:
0523           last = last - leng;
0524         }  // for( int iCl=LastCl
0525
0526       }  // if( nCl > 1
0527
0528     }  // if( (ia-vhit
0529
0530   new_ich:
0531     continue;
0532   }  //  for( ich=1
0533
0534   if (nCl > 0)
0535   {
0536     ib = 0;
0537     for (int iCl = 0; iCl < nCl; iCl++)
0538     {
0539       leng = LenCl[iCl];
0540       hl.clear();
0541       for (ich = 0; ich < leng; ich++)
0542       {
0543         hl.push_back(vhit[ib + ich]);// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
0544       }
0545       Clt.ReInitialize(hl);
0546       ib += LenCl[iCl];
0547       fClusters->push_back(Clt);
0548     }
0549   }
0550   delete[] LenCl;
0551   delete[] vhit;
0552   delete[] vt;
0553
0554   return nCl;
0555 }
0556
0557 // ///////////////////////////////////////////////////////////////////////////
0558
0559 void BEmcRec::Momenta(std::vector<EmcModule>* phit, float& pe, float& px,
0560                       float& py, float& pxx, float& pyy, float& pyx,
0561                       float thresh) const
0562 {
0563   // First and second momenta calculation
0564
0565   float a;
0566   float x;
0567   float y;
0568   float e;
0569   float xx;
0570   float yy;
0571   float yx;
0572   std::vector<EmcModule>::iterator ph;
0573
0574   pe = 0;
0575   px = 0;
0576   py = 0;
0577   pxx = 0;
0578   pyy = 0;
0579   pyx = 0;
0580   if (phit->empty())
0581   {
0582     return;
0583   }
0584
0585   // Find max energy tower
0586   //
0587   ph = phit->begin();
0588   float emax = 0;
0589   int ichmax = 0;
0590
0591   while (ph != phit->end())
0592   {
0593     a = ph->amp;
0594     if (a > emax)
0595     {
0596       emax = a;
0597       ichmax = ph->ich;
0598     }
0599     ++ph;
0600   }
0601
0602   if (emax <= 0)
0603   {
0604     return;
0605   }
0606
0607   int iymax = ichmax / fNx;
0608   int ixmax = ichmax - iymax * fNx;
0609
0610   // Calculate CG relative to max energy tower
0611
0612   x = 0;
0613   y = 0;
0614   e = 0;
0615   xx = 0;
0616   yy = 0;
0617   yx = 0;
0618   ph = phit->begin();
0619
0620   while (ph != phit->end())
0621   {
0622     a = ph->amp;
0623     if (a > thresh)
0624     {
0625       int iy = ph->ich / fNx;
0626       int ix = ph->ich - (iy * fNx);
0627       int idx = iTowerDist(ixmax, ix);
0628       int idy = iy - iymax;
0629       e += a;
0630       x += idx * a;
0631       y += idy * a;
0632       xx += a * idx * idx;
0633       yy += a * idy * idy;
0634       yx += a * idx * idy;
0635     }
0636     ++ph;
0637   }
0638
0639   pe = e;
0640
0641   if (e > 0)
0642   {
0643     x /= e;
0644     y /= e;
0645     xx = xx / e - x * x;
0646     yy = yy / e - y * y;
0647     yx = yx / e - y * x;
0648
0649     x += ixmax;
0650     y += iymax;
0651
0652     while (x < -0.5)
0653     {
0654       x += float(fNx);
0655     }
0656     while (x >= fNx - 0.5)
0657     {
0658       x -= float(fNx);
0659     }
0660
0661     px = x;
0662     py = y;
0663     pxx = xx;
0664     pyy = yy;
0665     pyx = yx;
0666   }
0667 }
0668
0669 // ///////////////////////////////////////////////////////////////////////////
0670
0671 float BEmcRec::PredictEnergy(float en, float xcg, float ycg, int ix, int iy)
0672 {
0673   if (_emcprof != nullptr && bProfileProb)
0674   {
0675     return PredictEnergyProb(en, xcg, ycg, ix, iy);
0676   }
0677
0678   float dx = std::fabs(fTowerDist(float(ix), xcg));
0679   float dy = ycg - iy;
0680   return PredictEnergyParam(en, dx, dy);
0681 }
0682
0683 float BEmcRec::PredictEnergyParam(float /*en*/, float xc, float yc)
0684 {
0685   // Calculates the energy deposited in the tower, the distance between
0686   // its center and shower Center of Gravity being (xc,yc)
0687   // en - shower energy
0688
0689   float dx;
0690   float dy;
0691   float r1;
0692   float r2;
0693   float r3;
0694   float fPpar1;
0695   float fPpar2;
0696   float fPpar3;
0697   float fPpar4;
0698
0699   float fPshiftx = 0;  // !!!!! Untill tuned ... may not be necessary
0700   float fPshifty = 0;  // !!!!! Untill tuned ... may not be necessary
0701
0702   /*
0703   float lgE;
0704   if (en <= 1.e-10)
0705     lgE = 0;
0706   else
0707     lgE = log(en);
0708   fPpar1=0.59-(1.45+0.13*lgE)*sin2a;
0709   fPpar2=0.265+(0.80+0.32*lgE)*sin2a;
0710   fPpar3=0.25+(0.45-0.036*lgE)*sin2a;
0711   fPpar4=0.42;
0712   */
0713   fPpar1 = 0.549;
0714   fPpar2 = 0.304;
0715   fPpar3 = 0.389;
0716   fPpar4 = 0.326;
0717   /*
0718   fPpar1 = 0.486;
0719   fPpar2 = 0.302;
0720   fPpar3 = 0.354;
0721   fPpar4 = 0.407;
0722   */
0723   /*
0724   fPpar1 = 0.343;
0725   fPpar2 = 0.509;
0726   fPpar3 = 0.199;
0727   fPpar4 = 0.548;
0728   */
0729
0730   //  if (en > 0) SetProfileParameters(-1, en, xc, yc);
0731
0732   dx = std::fabs(xc - fPshiftx);
0733   dy = std::fabs(yc - fPshifty);
0734   r2 = dx * dx + dy * dy;
0735   r1 = std::sqrt(r2);
0736   r3 = r2 * r1;
0737   double e = (fPpar1 * std::exp(-r3 / fPpar2)) + (fPpar3 * std::exp(-r1 / fPpar4));
0738
0739   return e;
0740 }
0741
0742 float BEmcRec::PredictEnergyProb(float en, float xcg, float ycg, int ix, int iy)
0743 // Predict tower energy from profiles used in GetProb()
0744 // This is expected to be used in BEmcCluster::GetSubClusters
0745 {
0746   if (_emcprof == nullptr)
0747   {
0748     return -1;
0749   }
0750
0751   while (xcg < -0.5)
0752   {
0753     xcg += float(fNx);
0754   }
0755   while (xcg >= fNx - 0.5)
0756   {
0757     xcg -= float(fNx);
0758   }
0759
0760 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
0761   int ixcg = int(xcg + 0.5);
0762 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
0763   int iycg = int(ycg + 0.5);
0764   float ddx = std::fabs(xcg - ixcg);
0765   float ddy = std::fabs(ycg - iycg);
0766
0767   float xg=0;
0768   float yg=0;
0769   float zg=0;
0770   Tower2Global(en, xcg, ycg, xg, yg, zg);
0771
0772   float theta;
0773   float phi;
0774   GetImpactThetaPhi(xg, yg, zg, theta, phi);
0775
0776   int isx = 1;
0777   if (xcg - ixcg < 0)
0778   {
0779     isx = -1;
0780   }
0781   int isy = 1;
0782   if (ycg - iycg < 0)
0783   {
0784     isy = -1;
0785   }
0786
0787   int idx = iTowerDist(ixcg, ix) * isx;
0788   int idy = (iy - iycg) * isy;
0789
0790   int id = -1;
0791   if (idx == 0 && idy == 0)
0792   {
0793     id = 0;
0794   }
0795   else if (idx == 1 && idy == 0)
0796   {
0797     id = 1;
0798   }
0799   else if (idx == 1 && idy == 1)
0800   {
0801     id = 2;
0802   }
0803   else if (idx == 0 && idy == 1)
0804   {
0805     id = 3;
0806   }
0807
0808   if (id < 0)
0809   {
0810     float dx = std::fabs(fTowerDist(xcg, float(ix)));
0811     float dy = std::fabs(iy - ycg);
0812     float rr = std::sqrt((dx * dx) + (dy * dy));
0813     //    return PredictEnergyParam(en, dx, dy);
0814     return _emcprof->PredictEnergyR(en, theta, phi, rr);
0815   }
0816
0817   float ep[4];
0818   float err[4];
0819   for (int ip = 0; ip < 4; ip++)
0820   {
0821     _emcprof->PredictEnergy(ip, en, theta, phi, ddx, ddy, ep[ip], err[ip]);
0822   }
0823
0824   float eout;
0825
0826   if (id == 0)
0827   {
0828     eout = (ep[1] + ep[2]) / 2. + ep[3];
0829   }
0830   else if (id == 1)
0831   {
0832     eout = (ep[0] - ep[2]) / 2. - ep[3];
0833   }
0834   else if (id == 3)
0835   {
0836     eout = (ep[0] - ep[1]) / 2. - ep[3];
0837   }
0838   else
0839   {
0840     eout = ep[3];
0841   }
0842
0843   //  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);
0844   if (eout < 0)
0845   {
0846     eout = 1e-6;
0847   }
0848
0849   return eout;
0850 }
0851
0852 // ///////////////////////////////////////////////////////////////////////////
0853
0854 float BEmcRec::GetTowerEnergy(int iy, int iz, std::vector<EmcModule>* plist) const
0855 {
0856   int nn = plist->size();
0857   if (nn <= 0)
0858   {
0859     return 0;
0860   }
0861
0862   for (int i = 0; i < nn; i++)
0863   {
0864     int ich = (*plist)[i].ich;
0865     int iyt = ich / fNx;
0866     int izt = ich % fNx;
0867     if (iy == iyt && iz == izt)
0868     {
0869       return (*plist)[i].amp;
0870     }
0871   }
0872   return 0;
0873 }
0874
0875 // !!!!! Change here to a ponter to HitList
0876 float BEmcRec::GetProb(std::vector<EmcModule> HitList, float en, float xg, float yg, float zg, float& chi2, int& ndf)
0877 // Do nothing; should be defined in a detector specific module BEmcRec{Name}
0878 {
0879   //  float enoise = 0.01;  // 10 MeV per tower
0880   float enoise = GetProbNoiseParam();
0881   //  float thresh = 0.01;
0882   float thresh = GetTowerThreshold();
0883
0884   chi2 = 0;
0885   ndf = 0;
0886   if (_emcprof == nullptr)
0887   {
0888     return -1;
0889   }
0890
0891   if (!(_emcprof->IsLoaded()))
0892   {
0893     return -1;
0894   }
0895
0896   int nn = HitList.size();
0897   if (nn <= 0)
0898   {
0899     return -1;
0900   }
0901
0902   float theta;
0903   float phi;
0904   GetImpactThetaPhi(xg, yg, zg, theta, phi);
0905
0906   // z coordinate below means x coordinate
0907
0908   float etot;
0909   float zcg;
0910   float ycg;
0911   float zz;
0912   float yy;
0913   float yz;
0914   Momenta(&HitList, etot, zcg, ycg, zz, yy, yz, thresh);
0915
0916 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
0917   int iz0cg = int(zcg + 0.5);
0918 // NOLINTNEXTLINE(bugprone-incorrect-roundings)
0919   int iy0cg = int(ycg + 0.5);
0920   float ddz = std::fabs(zcg - iz0cg);
0921   float ddy = std::fabs(ycg - iy0cg);
0922
0923   int isz = 1;
0924   if (zcg - iz0cg < 0)
0925   {
0926     isz = -1;
0927   }
0928   int isy = 1;
0929   if (ycg - iy0cg < 0)
0930   {
0931     isy = -1;
0932   }
0933
0934   // 4 central towers: 43
0935   //                   12
0936   // Tower 1 - central one
0937   float e1;
0938   float e2;
0939   float e3;
0940   float e4;
0941   e1 = GetTowerEnergy(iy0cg, iz0cg, &HitList);
0942   e2 = GetTowerEnergy(iy0cg, iz0cg + isz, &HitList);
0943   e3 = GetTowerEnergy(iy0cg + isy, iz0cg + isz, &HitList);
0944   e4 = GetTowerEnergy(iy0cg + isy, iz0cg, &HitList);
0945
0946   if (e1 < thresh)
0947   {
0948     e1 = 0;
0949   }
0950   if (e2 < thresh)
0951   {
0952     e2 = 0;
0953   }
0954   if (e3 < thresh)
0955   {
0956     e3 = 0;
0957   }
0958   if (e4 < thresh)
0959   {
0960     e4 = 0;
0961   }
0962
0963   float e1t = (e1 + e2 + e3 + e4) / etot;
0964   float e2t = (e1 + e2 - e3 - e4) / etot;
0965   float e3t = (e1 - e2 - e3 + e4) / etot;
0966   float e4t = (e3) / etot;
0967   //  float rr = sqrt((0.5-ddz)*(0.5-ddz)+(0.5-ddy)*(0.5-ddy));
0968
0969   // Predicted values
0970   const int NP = 4;  // From BEmcProfile
0971   float ep[NP];
0972   float err[NP];
0973   for (int ip = 0; ip < NP; ip++)
0974   {
0975     _emcprof->PredictEnergy(ip, en, theta, phi, ddz, ddy, ep[ip], err[ip]);
0976     if (ep[ip] < 0)
0977     {
0978       return -1;
0979     }
0980     if (ip < 3)
0981     {
0982       err[ip] = std::sqrt((err[ip] * err[ip]) + (4 * enoise * enoise / etot / etot));
0983     }
0984     else
0985     {
0986       err[ip] = std::sqrt((err[ip] * err[ip]) + (1 * enoise * enoise / etot / etot));
0987     }
0988   }
0989
0990   chi2 = 0.;
0991   chi2 += (ep[0] - e1t) * (ep[0] - e1t) / err[0] / err[0];
0992   chi2 += (ep[1] - e2t) * (ep[1] - e2t) / err[1] / err[1];
0993   chi2 += (ep[2] - e3t) * (ep[2] - e3t) / err[2] / err[2];
0994   chi2 += (ep[3] - e4t) * (ep[3] - e4t) / err[3] / err[3];
0995   ndf = 4;
0996
0997   chi2 /= 1.5;
0998
0999   float prob = TMath::Prob(chi2, ndf);
1000
1001   return prob;
1002 }
1003
1004 // ///////////////////////////////////////////////////////////////////////////
1005 // Static functions
1006
1007 int BEmcRec::HitNCompare(const void* h1, const void* h2)
1008 {
1009   return (static_cast<const EmcModule*>(h1)->ich - static_cast<const EmcModule*>(h2)->ich);
1010 }
1011
1012 // ///////////////////////////////////////////////////////////////////////////
1013
1014 int BEmcRec::HitACompare(const void* h1, const void* h2)
1015 {
1016   float amp1 = static_cast<const EmcModule*>(h1)->amp;
1017   float amp2 = static_cast<const EmcModule*>(h2)->amp;
1018   return (amp1 < amp2) ? 1 : (amp1 > amp2) ? -1 : 0; // NOLINT(readability-avoid-nested-conditional-operator)
1019 }
1020
1021 // ///////////////////////////////////////////////////////////////////////////
1022
1023 void BEmcRec::ZeroVector(int* v, int N)
1024 {
1025   int* p = v;
1026   for (int i = 0; i < N; i++)
1027   {
1028     *p++ = 0;
1029   }
1030 }
1031
1032 // ///////////////////////////////////////////////////////////////////////////
1033
1034 void BEmcRec::ZeroVector(float* v, int N)
1035 {
1036   float* p = v;
1037   for (int i = 0; i < N; i++)
1038   {
1039     *p++ = 0;
1040   }
1041 }
1042
1043 // ///////////////////////////////////////////////////////////////////////////
1044
1045 void BEmcRec::ZeroVector(EmcModule* v, int N)
1046 {
1047   //  memset(v, 0, N*sizeof(EmcModule));
1048   for (int i = 0; i < N; i++)
1049   {
1050     v[i].ich = 0;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
1051     v[i].amp = 0;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
1052     v[i].tof = 0;// NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
1053   }
1054 }
1055
1056 // ///////////////////////////////////////////////////////////////////////////
1057
1058 void BEmcRec::CopyVector(const int* from, int* to, int N)
1059 {
1060   if (N <= 0)
1061   {
1062     return;
1063   }
1064   for (int i = 0; i < N; i++)
1065   {
1066     to[i] = from[i];
1067   }
1068 }
1069
1070 // ///////////////////////////////////////////////////////////////////////////
1071
1072 void BEmcRec::CopyVector(const EmcModule* from, EmcModule* to, int N)
1073 {
1074   if (N <= 0)
1075   {
1076     return;
1077   }
1078   for (int i = 0; i < N; i++)
1079   {
1080     to[i] = from[i]; // NOLINT(bugprone-pointer-arithmetic-on-polymorphic-object)
1081   }
1082 }
1083
1084 // ///////////////////////////////////////////////////////////////////////////
1085
1086 /* Future improvements:
1087
1088 1. FindClusters(): to ensure that all EmcModules are above energy threshold
1089 set by SetThreshold routine (or default one)
1090
1091 */
1092
1093 // ///////////////////////////////////////////////////////////////////////////
1094 // EOF