Execute on an entire collection of clusters.
289 typedef std::pair<clusterIdx_t, clusterIdx_t> clusterPair_t;
290 std::vector<clusterPair_t> clusterIdx;
312 if (theClusColl->size() >= noCluster) {
313 msg(MSG::ERROR) <<
"Too many clusters" <<
endmsg;
314 return StatusCode::FAILURE;
319 clusterPair_t(noCluster, noCluster));
326 for(; cellIter != cellIterEnd; cellIter++ ){
333 if(hashid >= (signalCells)->
size())
continue;
334 myCell = (*signalCells).findCell(hashid);
335 if(!myCell)
continue;
340 if ( clusterIdx[(
unsigned int)myHashId].
first != noCluster) {
344 clusterIdx[(
unsigned int)myHashId].
first = iCluster;
347 clusterIdx[(
unsigned int)myHashId].
first = iCluster;
358 std::vector<CaloClusterMomentsMaker_DigiHSTruth_detail::cellinfo> cellinfo;
362 std::vector<IdentifierHash> theNeighbors;
367 for( ;clusIter!=clusIterEnd;++clusIter,++iClus) {
370 double w(0),xc(0),yc(0),zc(0);
371 double eBad(0),ebad_dac(0),ePos(0),eBadLArQ(0),sumSig2(0),maxAbsSig(0);
372 double eLAr2(0),eLAr2Q(0);
373 double eTile2(0),eTile2Q(0);
375 int nbad(0),nbad_dac(0),nBadLArHV(0);
376 unsigned int ncell(0),
i,nSigSampl(0);
377 unsigned int theNumOfCells = theCluster->
size();
378 double theClusterEnergy = 0;
379 double theClusterAbsEnergy = 0;
380 double theClusterEta = 0;
381 double theClusterPhi = 0;
385 int iCellScndMax(-1);
387 if (cellinfo.capacity() == 0)
388 cellinfo.reserve (theNumOfCells*2);
389 cellinfo.resize (theNumOfCells);
391 double phi0 = theCluster->
phi();;
396 std::fill (myMoments.begin(), myMoments.end(), 0);
397 std::fill (myNorms.begin(), myNorms.end(), 0);
406 for(; cellIter != cellIterEnd; cellIter++ ){
408 const CaloCell* pCell = (*cellIter);
413 if(hashid >= (signalCells)->
size())
continue;
414 pCell = (*signalCells).findCell(hashid);
419 double ene = pCell->
e();
424 double cellPhi = myCDDE->
phi();
427 theClusterEnergy +=
weight * ene;
428 theClusterAbsEnergy +=
weight*std::abs(ene);
429 theClusterEta +=
weight*std::abs(ene)*pCell->
eta();
430 theClusterPhi +=
weight*std::abs(ene)* thePhi;
442 && !((pCell->
provenance() & 0x0800) == 0x0800)) {
454 if ( ((tq1&0xFF) != 0xFF) && ((tq2&0xFF) != 0xFF) ) {
467 noise->getEffectiveSigma(pCell->
ID(),pCell->
gain(),pCell->
energy()) : \
472 if ( std::abs(Sig) > std::abs(maxAbsSig) ) {
482 if ( clusterIdx[myHashId].
first == iClus ) {
483 theNeighbors.clear();
485 for (
const auto& nhash: theNeighbors) {
486 clusterPair_t&
idx = clusterIdx[nhash];
489 if (
idx.second == iClus )
continue;
492 if (
idx.first == noCluster ) {
494 }
else if (
idx.first != iClus ) {
502 if ( myCDDE && ene > 0. &&
weight > 0) {
516 if (iCellMax < 0 || ci.energy > cellinfo[iCellMax].
energy ) {
517 iCellScndMax = iCellMax;
520 else if (iCellScndMax < 0 ||
521 ci.
energy > cellinfo[iCellScndMax].energy )
523 iCellScndMax =
ncell;
563 C(0,0) +=
e2*(ci.
x-xc)*(ci.
x-xc);
564 C(1,0) +=
e2*(ci.
x-xc)*(ci.
y-yc);
565 C(2,0) +=
e2*(ci.
x-xc)*(ci.
z-zc);
567 C(1,1) +=
e2*(ci.
y-yc)*(ci.
y-yc);
568 C(2,1) +=
e2*(ci.
y-yc)*(ci.
z-zc);
570 C(2,2) +=
e2*(ci.
z-zc)*(ci.
z-zc);
577 Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigensolver(
C);
578 if (eigensolver.info() != Eigen::Success) {
579 msg(MSG::WARNING) <<
"Failed to compute Eigenvalues -> Can't determine shower axis" <<
endmsg;
586 const Eigen::Vector3d&
S=eigensolver.eigenvalues();
587 const Eigen::Matrix3d& U=eigensolver.eigenvectors();
589 const double epsilon = 1.E-6;
591 if ( std::abs(
S[0]) >= epsilon && std::abs(
S[1]) >= epsilon && std::abs(
S[2]) >= epsilon ) {
601 double tmpAngle=
Amg::angle(tmpAxis,showerAxis);
603 if ( tmpAngle > 90*
deg ) {
604 tmpAngle = 180*
deg - tmpAngle;
608 if ( iEigen == -1 || tmpAngle <
angle ) {
620 deltaTheta = showerAxis.theta() - prAxis.theta();
629 << prAxis[
Amg::y] <<
", " << prAxis[
Amg::z] <<
") deviates more than "
631 <<
" deg from IP-to-ClusterCenter-axis (" << showerAxis[
Amg::x] <<
", "
632 << showerAxis[
Amg::y] <<
", " << showerAxis[
Amg::z] <<
")");
638 << showerAxis.y() <<
", " << showerAxis.z() <<
")");
644 for (
auto& ci : cellinfo) {
647 ci.
r = ((currentCell-showerCenter).cross(showerAxis)).
mag();
649 ci.
lambda = (currentCell-showerCenter).
dot(showerAxis);
655 double commonNorm = 0;
656 double phi0 =
ncell > 0 ? cellinfo[0].phi : 0;
667 case FIRST_ETA_DigiHSTruth:
670 case FIRST_PHI_DigiHSTruth:
676 case SECOND_R_DigiHSTruth:
677 myMoments[iMoment] += ci.
energy*ci.
r*ci.
r;
679 case SECOND_LAMBDA_DigiHSTruth:
682 case LATERAL_DigiHSTruth:
683 if ( (
int)
i != iCellMax && (
int)
i != iCellScndMax ) {
684 myMoments[iMoment] += ci.
energy*ci.
r*ci.
r;
685 myNorms[iMoment] += ci.
energy*ci.
r*ci.
r;
691 myNorms[iMoment] += rm*rm*ci.
energy;
694 case LONGITUDINAL_DigiHSTruth:
695 if ( (
int)
i != iCellMax && (
int)
i != iCellScndMax ) {
703 myNorms[iMoment] += lm*lm*ci.
energy;
706 case FIRST_ENG_DENS_DigiHSTruth:
709 myNorms[iMoment] += ci.
energy;
712 case SECOND_ENG_DENS_DigiHSTruth:
715 myNorms[iMoment] += ci.
energy;
718 case ENG_FRAC_EM_DigiHSTruth:
726 myMoments[iMoment] += ci.
energy;
728 case ENG_FRAC_MAX_DigiHSTruth:
729 if ( (
int)
i == iCellMax )
730 myMoments[iMoment] = ci.
energy;
740 eBadLArHV= hvFrac.first;
741 nBadLArHV=hvFrac.second;
751 case FIRST_ETA_DigiHSTruth:
752 case FIRST_PHI_DigiHSTruth:
753 case SECOND_R_DigiHSTruth:
754 case SECOND_LAMBDA_DigiHSTruth:
755 case ENG_FRAC_EM_DigiHSTruth:
756 case ENG_FRAC_MAX_DigiHSTruth:
757 myNorms[iMoment] = commonNorm;
759 case DELTA_PHI_DigiHSTruth:
762 case DELTA_THETA_DigiHSTruth:
763 myMoments[iMoment] = deltaTheta;
765 case DELTA_ALPHA_DigiHSTruth:
766 myMoments[iMoment] =
angle;
768 case CENTER_X_DigiHSTruth:
769 myMoments[iMoment] = showerCenter.x();
771 case CENTER_Y_DigiHSTruth:
772 myMoments[iMoment] = showerCenter.y();
774 case CENTER_Z_DigiHSTruth:
775 myMoments[iMoment] = showerCenter.z();
777 case CENTER_MAG_DigiHSTruth:
778 myMoments[iMoment] = showerCenter.mag();
780 case CENTER_LAMBDA_DigiHSTruth:
788 double r_calo(0),z_calo(0),lambda_c(0);
813 if ( z_calo != 0 && showerAxis.z() != 0 ) {
814 lambda_c = std::abs((z_calo-showerCenter.z())/showerAxis.z());
818 double r_s2 = showerAxis.x()*showerAxis.x()
819 +showerAxis.y()*showerAxis.y();
820 double r_cs = showerAxis.x()*showerCenter.x()
821 +showerAxis.y()*showerCenter.y();
822 double r_cr = showerCenter.x()*showerCenter.x()
823 +showerCenter.y()*showerCenter.y()-r_calo*r_calo;
825 double det = r_cs*r_cs/(r_s2*r_s2) - r_cr/r_s2;
828 double l1(-r_cs/r_s2);
832 if ( std::abs(
l1) < std::abs(
l2) )
833 lambda_c = std::abs(
l1);
835 lambda_c = std::abs(
l2);
839 myMoments[iMoment] = lambda_c;
842 case ENG_FRAC_CORE_DigiHSTruth:
844 myMoments[iMoment] += maxSampE[
i];
845 myNorms[iMoment] = commonNorm;
847 case ISOLATION_DigiHSTruth:
854 const double eSample = theCluster->
eSample(
s);
856 int nAll = nbEmpty[
i]+nbNonEmpty[
i];
858 myMoments[iMoment] += (eSample*nbEmpty[
i])/nAll;
859 myNorms[iMoment] += eSample;
866 case ENG_BAD_CELLS_DigiHSTruth:
867 myMoments[iMoment] = eBad;
869 case N_BAD_CELLS_DigiHSTruth:
870 myMoments[iMoment] = nbad;
872 case N_BAD_CELLS_CORR_DigiHSTruth:
873 myMoments[iMoment] = nbad_dac;
875 case BAD_CELLS_CORR_E_DigiHSTruth:
876 myMoments[iMoment] = ebad_dac;
878 case BADLARQ_FRAC_DigiHSTruth:
879 myMoments[iMoment] = eBadLArQ/(theCluster->
e()!=0.?theCluster->
e():1.);
881 case ENG_POS_DigiHSTruth:
882 myMoments[iMoment] = ePos;
884 case SIGNIFICANCE_DigiHSTruth:
885 myMoments[iMoment] = (sumSig2>0?theCluster->
e()/sqrt(sumSig2):0.);
887 case CELL_SIGNIFICANCE_DigiHSTruth:
888 myMoments[iMoment] = maxAbsSig;
890 case CELL_SIG_SAMPLING_DigiHSTruth:
891 myMoments[iMoment] = nSigSampl;
893 case AVG_LAR_Q_DigiHSTruth:
894 myMoments[iMoment] = eLAr2Q/(eLAr2>0?eLAr2:1);
896 case AVG_TILE_Q_DigiHSTruth:
897 myMoments[iMoment] = eTile2Q/(eTile2>0?eTile2:1);
899 case ENG_BAD_HV_CELLS_DigiHSTruth:
900 myMoments[iMoment] = eBadLArHV;
902 case N_BAD_HV_CELLS_DigiHSTruth:
903 myMoments[iMoment] = nBadLArHV;
905 case ENERGY_DigiHSTruth:
906 myMoments[iMoment] = theClusterEnergy;
908 case ETA_DigiHSTruth:
909 if(theClusterAbsEnergy > 0)
910 myMoments[iMoment] = theClusterEta / theClusterAbsEnergy;
912 myMoments[iMoment] = 0;
915 case PHI_DigiHSTruth:
916 if(theClusterAbsEnergy > 0)
919 myMoments[iMoment] = 0;
931 for (
size_t iMoment = 0; iMoment !=
size; ++iMoment) {
933 if ( myNorms[iMoment] != 0 )
934 myMoments[iMoment] /= myNorms[iMoment];
935 if ( moment == FIRST_PHI_DigiHSTruth )
943 return StatusCode::SUCCESS;