Execute on an entire collection of clusters.
344{
346
347
348
349 using clusterIdx_t = std::uint16_t;
350 typedef std::pair<clusterIdx_t, clusterIdx_t> clusterPair_t;
351 std::vector<clusterPair_t> clusterIdx;
352 const clusterIdx_t noCluster = std::numeric_limits<clusterIdx_t>::max();
353
354 const CaloNoise*
noise=
nullptr;
358 }
359
360 SG::ReadCondHandle<CaloDetDescrManager> caloMgrHandle{
m_caloMgrKey, ctx };
361 const CaloDetDescrManager* caloDDMgr = *caloMgrHandle;
362
363
364 int nbEmpty[CaloCell_ID::Unknown];
365 int nbNonEmpty[CaloCell_ID::Unknown];
366
367
368
369
371
372 if (theClusColl->size() >= noCluster) {
373 msg(MSG::ERROR) <<
"Too many clusters" <<
endmsg;
374 return StatusCode::FAILURE;
375 }
376
377
378 clusterIdx.resize(
m_calo_id->calo_cell_hash_max(),
379 clusterPair_t(noCluster, noCluster));
380
381 int iClus = 0;
383
386 for(; cellIter != cellIterEnd; cellIter++ ){
388 const CaloCell* pCell = *cellIter;
389
390 Identifier myId = pCell->
ID();
391 IdentifierHash myHashId =
m_calo_id->calo_cell_hash(myId);
392 if ( clusterIdx[(unsigned int)myHashId].first != noCluster) {
393
395 if ( weight > 0.5 )
396 clusterIdx[(
unsigned int)myHashId].first = iClus;
397 }
398 else {
399 clusterIdx[(
unsigned int)myHashId].first = iClus;
400 }
401 }
402 ++iClus;
403 }
404 }
405
406
407
408
409
410 std::vector<CaloClusterMomentsMaker_detail::cellinfo> cellinfo;
411 std::vector<double> maxSampE (CaloCell_ID::Unknown);
414 std::vector<std::tuple<int,int> > nCellsSamp; nCellsSamp.reserve(CaloCell_ID::Unknown);
415 std::vector<IdentifierHash> theNeighbors;
416
419 int iClus = 0;
420 for( ;clusIter!=clusIterEnd;++clusIter,++iClus) {
422
423 double w(0),xc(0),yc(0),zc(0),
mx(0),
my(0),
mz(0),
mass(0);
424 double eBad(0),ebad_dac(0),ePos(0),eBadLArQ(0),sumSig2(0),maxAbsSig(0);
425 double eLAr2(0),eLAr2Q(0);
426 double eTile2(0),eTile2Q(0);
427 double eBadLArHV(0);
428 int nbad(0),nbad_dac(0),nBadLArHV(0);
429 unsigned int i,nSigSampl(0);
430 unsigned int theNumOfCells = theCluster->
size();
431
432
433 int iCellMax(-1);
434 int iCellScndMax(-1);
435
436 cellinfo.clear();
437 if (cellinfo.capacity() == 0)
438 cellinfo.reserve (theNumOfCells*2);
439
440 for(i=0;
i<(
unsigned int)CaloCell_ID::Unknown;
i++)
441 maxSampE[i] = 0;
442
444 std::fill (myMoments.begin(), myMoments.end(), 0);
445 std::fill (myNorms.begin(), myNorms.end(), 0);
447 std::fill_n(nbNonEmpty, CaloCell_ID::Unknown, 0);
448 std::fill_n(nbEmpty, CaloCell_ID::Unknown, 0);
449 }
450
451
454 for(; cellIter != cellIterEnd; cellIter++ ){
456
457 const CaloCell* pCell = (*cellIter);
458 Identifier myId = pCell->
ID();
459 const CaloDetDescrElement* myCDDE = pCell->
caloDDE();
460 double ene = pCell->
e();
465 nbad++;
466 if(ene!=0){
468 nbad_dac++;
469 }
470 }
471 else {
472 if ( myCDDE && ! (myCDDE->
is_tile())
474 && !((pCell->
provenance() & 0x0800) == 0x0800)) {
477 }
480 }
481 if ( myCDDE && myCDDE->
is_tile() ) {
485
486 if ( ((tq1&0xFF) != 0xFF) && ((tq2&0xFF) != 0xFF) ) {
488
489
491 }
492 }
493 }
494 if ( ene > 0 ) {
496 }
497
499
501
505
507
511 if ( ( std::abs(Sig) > std::abs(maxAbsSig) ) ||
512 ( std::abs(Sig) == std::abs(maxAbsSig) && thisSampl > nSigSampl ) ||
513 ( std::abs(Sig) == std::abs(maxAbsSig) && thisSampl == nSigSampl && Sig > maxAbsSig ) ) {
514 maxAbsSig = Sig;
515 nSigSampl = thisSampl;
516 }
517
518 }
519 else {
520 if ( std::abs(Sig) > std::abs(maxAbsSig) ) {
521 maxAbsSig = Sig;
523 }
524 }
525 }
527
528
529
530 IdentifierHash myHashId =
m_calo_id->calo_cell_hash(myId);
531 if ( clusterIdx[myHashId].first == iClus ) {
532 theNeighbors.clear();
534 for (const auto& nhash: theNeighbors) {
535 clusterPair_t&
idx = clusterIdx[nhash];
536
537
538 if (
idx.second == iClus )
continue;
540
541 if (
idx.first == noCluster ) {
542 ++ nbEmpty[
m_calo_id->calo_sample(nhash)];
543 }
else if (
idx.first != iClus ) {
544 ++ nbNonEmpty[
m_calo_id->calo_sample(nhash)];
545 }
546
547 }
548 }
549 }
550
551 if ( myCDDE != nullptr ) {
554 size_t idx((
size_t)sam);
555 if ( idx >= nCellsSamp.size() ) { nCellsSamp.resize(idx+1, { 0, 0 } ); }
556 std::get<0>(nCellsSamp[idx])++;
557
558 if ( sam == CaloCell_ID::EME2 && std::abs(myCDDE->
eta()) >
m_etaInnerWheel ) { std::get<1>(nCellsSamp[idx])++; }
559 }
560 if ( ene > 0. && weight > 0) {
561
562 cellinfo.push_back (CaloClusterMomentsMaker_detail::cellinfo {
567 .eta = myCDDE->
eta(),
568 .phi = myCDDE->
phi(),
569 .r = 0,
570 .lambda = 0,
571 .volume = myCDDE->
volume(),
573 .identifier = cellIter.
index()
574
575
576
577 });
578
579 CaloClusterMomentsMaker_detail::cellinfo& ci = cellinfo.back();
580
583
585 if (iCellMax < 0 ||
586 ci.
energy > cellinfo[iCellMax].energy ||
587 (ci.
energy == cellinfo[iCellMax].energy && ci.
identifier > cellinfo[iCellMax].identifier) ) {
588 iCellScndMax = iCellMax;
589 iCellMax = cellinfo.size()-1;
590 }
591 else if (iCellScndMax < 0 ||
592 ci.
energy > cellinfo[iCellScndMax].energy ||
593 (ci.
energy == cellinfo[iCellScndMax].energy && ci.
identifier > cellinfo[iCellScndMax].identifier) )
594 {
595 iCellScndMax = cellinfo.size()-1;
596 }
597 }
598 else {
599 if (iCellMax < 0 || ci.energy > cellinfo[iCellMax].energy ) {
600 iCellScndMax = iCellMax;
601 iCellMax = cellinfo.size()-1;
602 }
603 else if (iCellScndMax < 0 ||
604 ci.
energy > cellinfo[iCellScndMax].energy )
605 {
606 iCellScndMax = cellinfo.size()-1;
607 }
608 }
609
613
614 double dir = ci.
x*ci.
x+ci.
y*ci.
y+ci.
z*ci.
z;
615
616 if ( dir > 0) {
619 }
623
625 }
626 }
627 }
630 eBadLArHV= hvFrac.first;
631 nBadLArHV=hvFrac.second;
632 }
633
636 if ( mass > 0) {
638 }
639 else {
640
642 }
643
644 if (w == 0) {
646 }
647
653
654
655
656
657
658
659
662
663
664
665
666
667
669 if ( cellinfo.size() > 2 ) {
670 Eigen::Matrix3d
C=Eigen::Matrix3d::Zero();
671 for(const CaloClusterMomentsMaker_detail::cellinfo& ci : cellinfo) {
673
674 C(0,0) +=
e2*(ci.
x-xc)*(ci.
x-xc);
675 C(1,0) +=
e2*(ci.
x-xc)*(ci.
y-yc);
676 C(2,0) +=
e2*(ci.
x-xc)*(ci.
z-zc);
677
678 C(1,1) +=
e2*(ci.
y-yc)*(ci.
y-yc);
679 C(2,1) +=
e2*(ci.
y-yc)*(ci.
z-zc);
680
681 C(2,2) +=
e2*(ci.
z-zc)*(ci.
z-zc);
683 }
684 C/=(
w != 0 ?
w : 1.0);
685
686 Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigensolver(
C);
687 if (eigensolver.info() != Eigen::Success) {
688 msg(MSG::WARNING) <<
"Failed to compute Eigenvalues -> Can't determine shower axis" <<
endmsg;
689 }
690 else {
691
692
693
694 const Eigen::Vector3d&
S=eigensolver.eigenvalues();
695 const Eigen::Matrix3d& U=eigensolver.eigenvectors();
696
698
699 if ( std::abs(S[0]) >= epsilon && std::abs(S[1]) >= epsilon && std::abs(S[2]) >= epsilon ) {
700
702 int iEigen = -1;
703
706
707
708 double tmpAngle=
Amg::angle(tmpAxis,showerAxis);
709
710 if ( tmpAngle > 90*
deg ) {
711 tmpAngle = 180*
deg - tmpAngle;
712 tmpAxis = -tmpAxis;
713 }
714
715 if ( iEigen == -1 || tmpAngle <
angle ) {
718 prAxis = tmpAxis;
719 }
720 }
721
722
723
725
726 deltaTheta = showerAxis.theta() - prAxis.theta();
727
728
729
731 showerAxis = prAxis;
732 }
733 else
735 << prAxis[
Amg::y] <<
", " << prAxis[
Amg::z] <<
") deviates more than "
737 <<
" deg from IP-to-ClusterCenter-axis (" << showerAxis[
Amg::x] <<
", "
738 << showerAxis[
Amg::y] <<
", " << showerAxis[
Amg::z] <<
")");
739 }
740 else {
741 ATH_MSG_DEBUG(
"Eigenvalues close to 0, do not use principal axis");
742 }
743 }
744 }
745
747 << showerAxis[
Amg::y] <<
", " << showerAxis[
Amg::z] <<
")");
748
749
750
751
752
753
754 for (auto& ci : cellinfo) {
756
757 ci.
r = ((currentCell-showerCenter).cross(showerAxis)).mag();
758
759 ci.
lambda = (currentCell-showerCenter).
dot(showerAxis);
760 }
761
762
763
764
765 double commonNorm = 0;
766 double phi0 = cellinfo.size() > 0 ? cellinfo[0].phi : 0;
767
768 for(
unsigned i=0;
i<cellinfo.size();
i++) {
769 const CaloClusterMomentsMaker_detail::cellinfo& ci = cellinfo[
i];
770
774 ++ iMoment)
775 {
776
778 case FIRST_ETA:
780 break;
781 case FIRST_PHI:
782
783
784
786 break;
788 myMoments[iMoment] += ci.
energy*ci.
r*ci.
r;
789 break;
792 break;
793 case LATERAL:
794 if ( (int)i != iCellMax && (int)i != iCellScndMax ) {
795 myMoments[iMoment] += ci.
energy*ci.
r*ci.
r;
796 myNorms[iMoment] += ci.
energy*ci.
r*ci.
r;
797 }
798 else {
802 myNorms[iMoment] += rm*rm*ci.
energy;
803 }
804 break;
805 case LONGITUDINAL:
806 if ( (int)i != iCellMax && (int)i != iCellScndMax ) {
809 }
810 else {
814 myNorms[iMoment] += lm*lm*ci.
energy;
815 }
816 break;
820 myNorms[iMoment] += ci.
energy;
821 }
822 break;
823 case SECOND_ENG_DENS:
826 myNorms[iMoment] += ci.
energy;
827 }
828 break;
829 case ENG_FRAC_EM:
830 if ( ci.
sample == CaloCell_ID::EMB1
831 || ci.
sample == CaloCell_ID::EMB2
832 || ci.
sample == CaloCell_ID::EMB3
833 || ci.
sample == CaloCell_ID::EME1
834 || ci.
sample == CaloCell_ID::EME2
835 || ci.
sample == CaloCell_ID::EME3
836 || ci.
sample == CaloCell_ID::FCAL0 )
837 myMoments[iMoment] += ci.
energy;
838 break;
839 case ENG_FRAC_MAX:
840 if ( (int)i == iCellMax )
841 myMoments[iMoment] = ci.
energy;
842 break;
843 case PTD:
844
845
846
848 myNorms[iMoment] += ci.
energy;
849 break;
850 default:
851
852 break;
853 }
854 }
855 }
856
857
858
861 ++ iMoment)
862 {
863
865 case FIRST_ETA:
866 case FIRST_PHI:
869 case ENG_FRAC_EM:
870 case ENG_FRAC_MAX:
871 myNorms[iMoment] = commonNorm;
872 break;
873 case DELTA_PHI:
875 break;
876 case DELTA_THETA:
877 myMoments[iMoment] = deltaTheta;
878 break;
879 case DELTA_ALPHA:
880 myMoments[iMoment] =
angle;
881 break;
882 case CENTER_X:
883 myMoments[iMoment] = showerCenter.x();
884 break;
885 case CENTER_Y:
886 myMoments[iMoment] = showerCenter.y();
887 break;
888 case CENTER_Z:
889 myMoments[iMoment] = showerCenter.z();
890 break;
892 myMoments[iMoment] = showerCenter.mag();
893 break;
895
896
897
898
899
900
901 {
902 double r_calo(0),z_calo(0),lambda_c(0);
904 showerCenter.eta(),
905 showerCenter.phi(),
906 caloDDMgr);
907 if ( r_calo == 0 ) {
909 showerCenter.eta(),
910 showerCenter.phi(),
911 caloDDMgr);
912 if ( z_calo == 0 )
914 showerCenter.eta(),
915 showerCenter.phi(),
916 caloDDMgr);
917 if ( z_calo == 0 )
919 showerCenter.eta(),
920 showerCenter.phi(),
921 caloDDMgr);
922 if ( z_calo == 0 )
924 showerCenter.eta(),
925 showerCenter.phi(),
926 caloDDMgr);
927 if ( z_calo != 0 && showerAxis.z() != 0 ) {
928 lambda_c = std::abs((z_calo-showerCenter.z())/showerAxis.z());
929 }
930 }
931 else {
932 double r_s2 = showerAxis.x()*showerAxis.x()
933 +showerAxis.y()*showerAxis.y();
934 double r_cs = showerAxis.x()*showerCenter.x()
935 +showerAxis.y()*showerCenter.y();
936 double r_cr = showerCenter.x()*showerCenter.x()
937 +showerCenter.y()*showerCenter.y()-r_calo*r_calo;
938 if ( r_s2 > 0 ) {
939 double det = r_cs*r_cs/(r_s2*r_s2) - r_cr/r_s2;
940 if ( det > 0 ) {
942 double l1(-r_cs/r_s2);
946 if ( std::abs(l1) < std::abs(l2) )
947 lambda_c = std::abs(l1);
948 else
949 lambda_c = std::abs(l2);
950 }
951 }
952 }
953 myMoments[iMoment] = lambda_c;
954 }
955 break;
956 case ENG_FRAC_CORE:
957 for(i=0;
i<(
int)CaloCell_ID::Unknown;
i++)
958 myMoments[iMoment] += maxSampE[i];
959 myNorms[iMoment] = commonNorm;
960 break;
961 case ISOLATION:
962 {
963
964
965 for(
unsigned int i=0;
i != CaloSampling::Unknown; ++
i) {
968 const double eSample = theCluster->
eSample(s);
969 if (eSample > 0) {
970 int nAll = nbEmpty[
i]+nbNonEmpty[
i];
971 if (nAll > 0) {
972 myMoments[iMoment] += (eSample*nbEmpty[
i])/nAll;
973 myNorms[iMoment] += eSample;
974 }
975 }
976 }
977 }
978 }
979 break;
980 case ENG_BAD_CELLS:
981 myMoments[iMoment] = eBad;
982 break;
983 case N_BAD_CELLS:
984 myMoments[iMoment] = nbad;
985 break;
986 case N_BAD_CELLS_CORR:
987 myMoments[iMoment] = nbad_dac;
988 break;
989 case BAD_CELLS_CORR_E:
990 myMoments[iMoment] = ebad_dac;
991 break;
992 case BADLARQ_FRAC:
993 myMoments[iMoment] = eBadLArQ/(theCluster->
e()!=0.?theCluster->
e():1.);
994 break;
995 case ENG_POS:
996 myMoments[iMoment] = ePos;
997 break;
998 case SIGNIFICANCE:
999 myMoments[iMoment] = (sumSig2>0?theCluster->
e()/sqrt(sumSig2):0.);
1000 break;
1001 case CELL_SIGNIFICANCE:
1002 myMoments[iMoment] = maxAbsSig;
1003 break;
1004 case CELL_SIG_SAMPLING:
1005 myMoments[iMoment] = nSigSampl;
1006 break;
1007 case AVG_LAR_Q:
1008 myMoments[iMoment] = eLAr2Q/(eLAr2>0?eLAr2:1);
1009 break;
1010 case AVG_TILE_Q:
1011 myMoments[iMoment] = eTile2Q/(eTile2>0?eTile2:1);
1012 break;
1013 case ENG_BAD_HV_CELLS:
1014 myMoments[iMoment] = eBadLArHV;
1015 break;
1016 case N_BAD_HV_CELLS:
1017 myMoments[iMoment] = nBadLArHV;
1018 break;
1019 case PTD:
1020 myMoments[iMoment] = sqrt(myMoments[iMoment]);
1021 break;
1022 case MASS:
1023 myMoments[iMoment] =
mass;
1024 break;
1025 default:
1026
1027 break;
1028 }
1029 }
1030 }
1031
1032
1034 for (
size_t iMoment = 0; iMoment !=
size; ++iMoment) {
1036 if ( myNorms[iMoment] != 0 )
1037 myMoments[iMoment] /= myNorms[iMoment];
1038 if ( moment == FIRST_PHI )
1041 }
1042 }
1043
1045
1047 for ( size_t isam(0); isam < nCellsSamp.size(); ++isam ) {
1049 if ( isam == (size_t)CaloCell_ID::EME2 && std::get<1>(nCellsSamp.at(isam)) > 0 ) {
1051 }
1052 }
1053 nCellsSamp.clear();
1054 }
1055 }
1056
1057 return StatusCode::SUCCESS;
1058}
Scalar deltaPhi(const MatrixBase< Derived > &vec) const
size_t size() const
Number of registered mappings.
double angle(const GeoTrf::Vector2D &a, const GeoTrf::Vector2D &b)
CaloSampling::CaloSample CaloSample
virtual double e() const override final
get energy (data member) (synonym to method energy()
double energy() const
get energy (data member)
const CaloDetDescrElement * caloDDE() const
get pointer to CaloDetDescrElement (data member)
uint16_t provenance() const
get provenance (data member)
uint16_t quality() const
get quality (data member)
virtual bool badcell() const
check is cell is dead
CaloGain::CaloGain gain() const
get gain (data member )
Identifier ID() const
get ID (from cached data member) non-virtual and inline for fast access
unsigned index() const
Accessor for the index of the cell in the CaloCellContainer.
weight_t weight() const
Accessor for weight associated to this cell.
SG::ReadCondHandleKey< CaloDetDescrManager > m_caloMgrKey
bool m_secondTime
Retrieve second moment of cell times and store as moment.
std::vector< xAOD::CaloCluster::MomentType > m_validMoments
set of moments which will be calculated.
Gaudi::Property< bool > m_useGPUCriteria
bool m_nCellsPerSampling
store number of cells per sampling layer as moment
SG::ReadCondHandleKey< CaloNoise > m_noiseCDOKey
Key of the CaloNoise Conditions data object.
bool is_tile() const
cell belongs to Tile
CaloCell_ID::CaloSample getSampling() const
cell sampling
float eta() const
cell eta
float phi() const
cell phi
float volume() const
cell volume
static double fix(double phi)
static double diff(double phi1, double phi2)
simple phi1 - phi2 calculation, but result is fixed to respect range.
DataModel_detail::iterator< DataVector > iterator
const CaloClusterCellLink * getCellLinks() const
Get a pointer to the CaloClusterCellLink object (const version).
size_t size() const
size method (forwarded from CaloClusterCellLink obj)
CaloClusterCellLink::iterator cell_iterator
Iterator of the underlying CaloClusterCellLink (non-const version).
virtual double e() const
The total energy of the particle.
void insertMoment(MomentType type, double value)
const_cell_iterator cell_end() const
float eSample(const CaloSample sampling) const
flt_t secondTime() const
Access second moment of cell timing distribution.
MomentType
Enums to identify different moments.
const_cell_iterator cell_begin() const
Iterator of the underlying CaloClusterCellLink (const version).
void setNumberCellsInSampling(CaloSampling::CaloSample samp, int ncells, bool isInnerWheel=false)
Set the number of cells in a sampling layer.
CaloSampling::CaloSample CaloSample
bool hasSampling(const CaloSample s) const
Checks if certain smapling contributes to cluster.
void setMag(Amg::Vector3D &v, double mag)
scales the vector length without changing the angles
double angle(const Amg::Vector3D &v1, const Amg::Vector3D &v2)
calculates the opening angle between two vectors
Eigen::Matrix< double, 3, 1 > Vector3D
void nextDDE(Iter iter, Iter endIter)
Prefetch next CaloDDE.
void prefetchNext(Iter iter, Iter endIter)
Prefetch next object in sequence.
bool FIRST_ENG_DENS(const xAOD::TauJet &, const xAOD::CaloVertexedTopoCluster &cluster, float &out)
bool CENTER_MAG(const xAOD::TauJet &, const xAOD::CaloVertexedTopoCluster &cluster, float &out)
bool CENTER_LAMBDA(const xAOD::TauJet &, const xAOD::CaloVertexedTopoCluster &cluster, float &out)
bool SECOND_R(const xAOD::TauJet &, const xAOD::CaloVertexedTopoCluster &cluster, float &out)
bool SECOND_LAMBDA(const xAOD::TauJet &, const xAOD::CaloVertexedTopoCluster &cluster, float &out)
dot(G, fn, nodesToHighlight=[])
double e2(const xAOD::CaloCluster &cluster)
return the uncorrected cluster energy in 2nd sampling
CaloCluster_v1 CaloCluster
Define the latest version of the calorimeter cluster class.
double proxim(double b, double a)
CaloCell_ID::CaloSample sample
#define CXXUTILS_TRAPPING_FP