Execute on an entire collection of clusters.
343{
345
346
347
348 using clusterIdx_t = std::uint16_t;
349 typedef std::pair<clusterIdx_t, clusterIdx_t> clusterPair_t;
350 std::vector<clusterPair_t> clusterIdx;
351 const clusterIdx_t noCluster = std::numeric_limits<clusterIdx_t>::max();
352
353 const CaloNoise*
noise=
nullptr;
357 }
358
359 SG::ReadCondHandle<CaloDetDescrManager> caloMgrHandle{
m_caloMgrKey, ctx };
360 const CaloDetDescrManager* caloDDMgr = *caloMgrHandle;
361
362
363 int nbEmpty[CaloCell_ID::Unknown];
364 int nbNonEmpty[CaloCell_ID::Unknown];
365
366
367
368
370
371 if (theClusColl->size() >= noCluster) {
372 msg(MSG::ERROR) <<
"Too many clusters" <<
endmsg;
373 return StatusCode::FAILURE;
374 }
375
376
377 clusterIdx.resize(
m_calo_id->calo_cell_hash_max(),
378 clusterPair_t(noCluster, noCluster));
379
380 int iClus = 0;
382
385 for(; cellIter != cellIterEnd; cellIter++ ){
387 const CaloCell* pCell = *cellIter;
388
389 Identifier myId = pCell->
ID();
390 IdentifierHash myHashId =
m_calo_id->calo_cell_hash(myId);
391 if ( clusterIdx[(unsigned int)myHashId].first != noCluster) {
392
394 if ( weight > 0.5 )
395 clusterIdx[(
unsigned int)myHashId].first = iClus;
396 }
397 else {
398 clusterIdx[(
unsigned int)myHashId].first = iClus;
399 }
400 }
401 ++iClus;
402 }
403 }
404
405
406
407
408
409 std::vector<CaloClusterMomentsMaker_detail::cellinfo> cellinfo;
410 std::vector<double> maxSampE (CaloCell_ID::Unknown);
413 std::vector<std::tuple<int,int> > nCellsSamp; nCellsSamp.reserve(CaloCell_ID::Unknown);
414 std::vector<IdentifierHash> theNeighbors;
415
418 int iClus = 0;
419 for( ;clusIter!=clusIterEnd;++clusIter,++iClus) {
421
422 double w(0),xc(0),yc(0),zc(0),
mx(0),
my(0),
mz(0),
mass(0);
423 double eBad(0),ebad_dac(0),ePos(0),eBadLArQ(0),sumSig2(0),maxAbsSig(0);
424 double eLAr2(0),eLAr2Q(0);
425 double eTile2(0),eTile2Q(0);
426 double eBadLArHV(0);
427 int nbad(0),nbad_dac(0),nBadLArHV(0);
428 unsigned int i,nSigSampl(0);
429 unsigned int theNumOfCells = theCluster->
size();
430
431
432 int iCellMax(-1);
433 int iCellScndMax(-1);
434
435 cellinfo.clear();
436 if (cellinfo.capacity() == 0)
437 cellinfo.reserve (theNumOfCells*2);
438
439 for(i=0;
i<(
unsigned int)CaloCell_ID::Unknown;
i++)
440 maxSampE[i] = 0;
441
443 std::fill (myMoments.begin(), myMoments.end(), 0);
444 std::fill (myNorms.begin(), myNorms.end(), 0);
446 std::fill_n(nbNonEmpty, CaloCell_ID::Unknown, 0);
447 std::fill_n(nbEmpty, CaloCell_ID::Unknown, 0);
448 }
449
450
453 for(; cellIter != cellIterEnd; cellIter++ ){
455
456 const CaloCell* pCell = (*cellIter);
457 Identifier myId = pCell->
ID();
458 const CaloDetDescrElement* myCDDE = pCell->
caloDDE();
459 double ene = pCell->
e();
464 nbad++;
465 if(ene!=0){
467 nbad_dac++;
468 }
469 }
470 else {
471 if ( myCDDE && ! (myCDDE->
is_tile())
473 && !((pCell->
provenance() & 0x0800) == 0x0800)) {
476 }
479 }
480 if ( myCDDE && myCDDE->
is_tile() ) {
484
485 if ( ((tq1&0xFF) != 0xFF) && ((tq2&0xFF) != 0xFF) ) {
487
488
490 }
491 }
492 }
493 if ( ene > 0 ) {
495 }
496
501
503
507 if ( ( std::abs(Sig) > std::abs(maxAbsSig) ) ||
508 ( std::abs(Sig) == std::abs(maxAbsSig) && thisSampl > nSigSampl ) ||
509 ( std::abs(Sig) == std::abs(maxAbsSig) && thisSampl == nSigSampl && Sig > maxAbsSig ) ) {
510 maxAbsSig = Sig;
511 nSigSampl = thisSampl;
512 }
513
514 }
515 else {
516 if ( std::abs(Sig) > std::abs(maxAbsSig) ) {
517 maxAbsSig = Sig;
519 }
520 }
521 }
523
524
525
526 IdentifierHash myHashId =
m_calo_id->calo_cell_hash(myId);
527 if ( clusterIdx[myHashId].first == iClus ) {
528 theNeighbors.clear();
530 for (const auto& nhash: theNeighbors) {
531 clusterPair_t&
idx = clusterIdx[nhash];
532
533
534 if (
idx.second == iClus )
continue;
536
537 if (
idx.first == noCluster ) {
538 ++ nbEmpty[
m_calo_id->calo_sample(nhash)];
539 }
else if (
idx.first != iClus ) {
540 ++ nbNonEmpty[
m_calo_id->calo_sample(nhash)];
541 }
542
543 }
544 }
545 }
546
547 if ( myCDDE != nullptr ) {
550 size_t idx((
size_t)sam);
551 if ( idx >= nCellsSamp.size() ) { nCellsSamp.resize(idx+1, { 0, 0 } ); }
552 std::get<0>(nCellsSamp[idx])++;
553
554 if ( sam == CaloCell_ID::EME2 && std::abs(myCDDE->
eta()) >
m_etaInnerWheel ) { std::get<1>(nCellsSamp[idx])++; }
555 }
556 if ( ene > 0. && weight > 0) {
557
558 cellinfo.push_back (CaloClusterMomentsMaker_detail::cellinfo {
563 .eta = myCDDE->
eta(),
564 .phi = myCDDE->
phi(),
565 .r = 0,
566 .lambda = 0,
567 .volume = myCDDE->
volume(),
569 .identifier = cellIter.
index()
570
571
572
573 });
574
575 CaloClusterMomentsMaker_detail::cellinfo& ci = cellinfo.back();
576
579
581 if (iCellMax < 0 ||
582 ci.
energy > cellinfo[iCellMax].energy ||
583 (ci.
energy == cellinfo[iCellMax].energy && ci.
identifier > cellinfo[iCellMax].identifier) ) {
584 iCellScndMax = iCellMax;
585 iCellMax = cellinfo.size()-1;
586 }
587 else if (iCellScndMax < 0 ||
588 ci.
energy > cellinfo[iCellScndMax].energy ||
589 (ci.
energy == cellinfo[iCellScndMax].energy && ci.
identifier > cellinfo[iCellScndMax].identifier) )
590 {
591 iCellScndMax = cellinfo.size()-1;
592 }
593 }
594 else {
595 if (iCellMax < 0 || ci.energy > cellinfo[iCellMax].energy ) {
596 iCellScndMax = iCellMax;
597 iCellMax = cellinfo.size()-1;
598 }
599 else if (iCellScndMax < 0 ||
600 ci.
energy > cellinfo[iCellScndMax].energy )
601 {
602 iCellScndMax = cellinfo.size()-1;
603 }
604 }
605
609
610 double dir = ci.
x*ci.
x+ci.
y*ci.
y+ci.
z*ci.
z;
611
612 if ( dir > 0) {
615 }
619
621 }
622 }
623 }
626 eBadLArHV= hvFrac.first;
627 nBadLArHV=hvFrac.second;
628 }
629
632 if ( mass > 0) {
634 }
635 else {
636
638 }
639
640 if (w == 0) {
642 }
643
649
650
651
652
653
654
655
658
659
660
661
662
663
665 if ( cellinfo.size() > 2 ) {
666 Eigen::Matrix3d
C=Eigen::Matrix3d::Zero();
667 for(const CaloClusterMomentsMaker_detail::cellinfo& ci : cellinfo) {
669
670 C(0,0) +=
e2*(ci.
x-xc)*(ci.
x-xc);
671 C(1,0) +=
e2*(ci.
x-xc)*(ci.
y-yc);
672 C(2,0) +=
e2*(ci.
x-xc)*(ci.
z-zc);
673
674 C(1,1) +=
e2*(ci.
y-yc)*(ci.
y-yc);
675 C(2,1) +=
e2*(ci.
y-yc)*(ci.
z-zc);
676
677 C(2,2) +=
e2*(ci.
z-zc)*(ci.
z-zc);
679 }
680 C/=(
w != 0 ?
w : 1.0);
681
682 Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigensolver(
C);
683 if (eigensolver.info() != Eigen::Success) {
684 msg(MSG::WARNING) <<
"Failed to compute Eigenvalues -> Can't determine shower axis" <<
endmsg;
685 }
686 else {
687
688
689
690 const Eigen::Vector3d&
S=eigensolver.eigenvalues();
691 const Eigen::Matrix3d& U=eigensolver.eigenvectors();
692
694
695 if ( std::abs(S[0]) >= epsilon && std::abs(S[1]) >= epsilon && std::abs(S[2]) >= epsilon ) {
696
698 int iEigen = -1;
699
702
703
704 double tmpAngle=
Amg::angle(tmpAxis,showerAxis);
705
706 if ( tmpAngle > 90*
deg ) {
707 tmpAngle = 180*
deg - tmpAngle;
708 tmpAxis = -tmpAxis;
709 }
710
711 if ( iEigen == -1 || tmpAngle <
angle ) {
714 prAxis = tmpAxis;
715 }
716 }
717
718
719
721
722 deltaTheta = showerAxis.theta() - prAxis.theta();
723
724
725
727 showerAxis = prAxis;
728 }
729 else
731 << prAxis[
Amg::y] <<
", " << prAxis[
Amg::z] <<
") deviates more than "
733 <<
" deg from IP-to-ClusterCenter-axis (" << showerAxis[
Amg::x] <<
", "
734 << showerAxis[
Amg::y] <<
", " << showerAxis[
Amg::z] <<
")");
735 }
736 else {
737 ATH_MSG_DEBUG(
"Eigenvalues close to 0, do not use principal axis");
738 }
739 }
740 }
741
743 << showerAxis[
Amg::y] <<
", " << showerAxis[
Amg::z] <<
")");
744
745
746
747
748
749
750 for (auto& ci : cellinfo) {
752
753 ci.
r = ((currentCell-showerCenter).cross(showerAxis)).mag();
754
755 ci.
lambda = (currentCell-showerCenter).
dot(showerAxis);
756 }
757
758
759
760
761 double commonNorm = 0;
762 double phi0 = cellinfo.size() > 0 ? cellinfo[0].phi : 0;
763
764 for(
unsigned i=0;
i<cellinfo.size();
i++) {
765 const CaloClusterMomentsMaker_detail::cellinfo& ci = cellinfo[
i];
766
769 iMoment != size;
770 ++ iMoment)
771 {
772
774 case FIRST_ETA:
776 break;
777 case FIRST_PHI:
778
779
780
782 break;
784 myMoments[iMoment] += ci.
energy*ci.
r*ci.
r;
785 break;
788 break;
789 case LATERAL:
790 if ( (int)i != iCellMax && (int)i != iCellScndMax ) {
791 myMoments[iMoment] += ci.
energy*ci.
r*ci.
r;
792 myNorms[iMoment] += ci.
energy*ci.
r*ci.
r;
793 }
794 else {
798 myNorms[iMoment] += rm*rm*ci.
energy;
799 }
800 break;
801 case LONGITUDINAL:
802 if ( (int)i != iCellMax && (int)i != iCellScndMax ) {
805 }
806 else {
810 myNorms[iMoment] += lm*lm*ci.
energy;
811 }
812 break;
816 myNorms[iMoment] += ci.
energy;
817 }
818 break;
819 case SECOND_ENG_DENS:
822 myNorms[iMoment] += ci.
energy;
823 }
824 break;
825 case ENG_FRAC_EM:
826 if ( ci.
sample == CaloCell_ID::EMB1
827 || ci.
sample == CaloCell_ID::EMB2
828 || ci.
sample == CaloCell_ID::EMB3
829 || ci.
sample == CaloCell_ID::EME1
830 || ci.
sample == CaloCell_ID::EME2
831 || ci.
sample == CaloCell_ID::EME3
832 || ci.
sample == CaloCell_ID::FCAL0 )
833 myMoments[iMoment] += ci.
energy;
834 break;
835 case ENG_FRAC_MAX:
836 if ( (int)i == iCellMax )
837 myMoments[iMoment] = ci.
energy;
838 break;
839 case PTD:
840
841
842
844 myNorms[iMoment] += ci.
energy;
845 break;
846 default:
847
848 break;
849 }
850 }
851 }
852
853
854
856 iMoment != size;
857 ++ iMoment)
858 {
859
861 case FIRST_ETA:
862 case FIRST_PHI:
865 case ENG_FRAC_EM:
866 case ENG_FRAC_MAX:
867 myNorms[iMoment] = commonNorm;
868 break;
869 case DELTA_PHI:
871 break;
872 case DELTA_THETA:
873 myMoments[iMoment] = deltaTheta;
874 break;
875 case DELTA_ALPHA:
876 myMoments[iMoment] =
angle;
877 break;
878 case CENTER_X:
879 myMoments[iMoment] = showerCenter.x();
880 break;
881 case CENTER_Y:
882 myMoments[iMoment] = showerCenter.y();
883 break;
884 case CENTER_Z:
885 myMoments[iMoment] = showerCenter.z();
886 break;
888 myMoments[iMoment] = showerCenter.mag();
889 break;
891
892
893
894
895
896
897 {
898 double r_calo(0),z_calo(0),lambda_c(0);
900 showerCenter.eta(),
901 showerCenter.phi(),
902 caloDDMgr);
903 if ( r_calo == 0 ) {
905 showerCenter.eta(),
906 showerCenter.phi(),
907 caloDDMgr);
908 if ( z_calo == 0 )
910 showerCenter.eta(),
911 showerCenter.phi(),
912 caloDDMgr);
913 if ( z_calo == 0 )
915 showerCenter.eta(),
916 showerCenter.phi(),
917 caloDDMgr);
918 if ( z_calo == 0 )
920 showerCenter.eta(),
921 showerCenter.phi(),
922 caloDDMgr);
923 if ( z_calo != 0 && showerAxis.z() != 0 ) {
924 lambda_c = std::abs((z_calo-showerCenter.z())/showerAxis.z());
925 }
926 }
927 else {
928 double r_s2 = showerAxis.x()*showerAxis.x()
929 +showerAxis.y()*showerAxis.y();
930 double r_cs = showerAxis.x()*showerCenter.x()
931 +showerAxis.y()*showerCenter.y();
932 double r_cr = showerCenter.x()*showerCenter.x()
933 +showerCenter.y()*showerCenter.y()-r_calo*r_calo;
934 if ( r_s2 > 0 ) {
935 double det = r_cs*r_cs/(r_s2*r_s2) - r_cr/r_s2;
936 if ( det > 0 ) {
938 double l1(-r_cs/r_s2);
942 if ( std::abs(l1) < std::abs(l2) )
943 lambda_c = std::abs(l1);
944 else
945 lambda_c = std::abs(l2);
946 }
947 }
948 }
949 myMoments[iMoment] = lambda_c;
950 }
951 break;
952 case ENG_FRAC_CORE:
953 for(i=0;
i<(
int)CaloCell_ID::Unknown;
i++)
954 myMoments[iMoment] += maxSampE[i];
955 myNorms[iMoment] = commonNorm;
956 break;
957 case ISOLATION:
958 {
959
960
961 for(
unsigned int i=0;
i != CaloSampling::Unknown; ++
i) {
964 const double eSample = theCluster->
eSample(s);
965 if (eSample > 0) {
966 int nAll = nbEmpty[
i]+nbNonEmpty[
i];
967 if (nAll > 0) {
968 myMoments[iMoment] += (eSample*nbEmpty[
i])/nAll;
969 myNorms[iMoment] += eSample;
970 }
971 }
972 }
973 }
974 }
975 break;
976 case ENG_BAD_CELLS:
977 myMoments[iMoment] = eBad;
978 break;
979 case N_BAD_CELLS:
980 myMoments[iMoment] = nbad;
981 break;
982 case N_BAD_CELLS_CORR:
983 myMoments[iMoment] = nbad_dac;
984 break;
985 case BAD_CELLS_CORR_E:
986 myMoments[iMoment] = ebad_dac;
987 break;
988 case BADLARQ_FRAC:
989 myMoments[iMoment] = eBadLArQ/(theCluster->
e()!=0.?theCluster->
e():1.);
990 break;
991 case ENG_POS:
992 myMoments[iMoment] = ePos;
993 break;
994 case SIGNIFICANCE:
995 myMoments[iMoment] = (sumSig2>0?theCluster->
e()/sqrt(sumSig2):0.);
996 break;
997 case CELL_SIGNIFICANCE:
998 myMoments[iMoment] = maxAbsSig;
999 break;
1000 case CELL_SIG_SAMPLING:
1001 myMoments[iMoment] = nSigSampl;
1002 break;
1003 case AVG_LAR_Q:
1004 myMoments[iMoment] = eLAr2Q/(eLAr2>0?eLAr2:1);
1005 break;
1006 case AVG_TILE_Q:
1007 myMoments[iMoment] = eTile2Q/(eTile2>0?eTile2:1);
1008 break;
1009 case ENG_BAD_HV_CELLS:
1010 myMoments[iMoment] = eBadLArHV;
1011 break;
1012 case N_BAD_HV_CELLS:
1013 myMoments[iMoment] = nBadLArHV;
1014 break;
1015 case PTD:
1016 myMoments[iMoment] = sqrt(myMoments[iMoment]);
1017 break;
1018 case MASS:
1019 myMoments[iMoment] =
mass;
1020 break;
1021 default:
1022
1023 break;
1024 }
1025 }
1026 }
1027
1028
1030 for (
size_t iMoment = 0; iMoment !=
size; ++iMoment) {
1032 if ( myNorms[iMoment] != 0 )
1033 myMoments[iMoment] /= myNorms[iMoment];
1034 if ( moment == FIRST_PHI )
1037 }
1038 }
1039
1041
1043 for ( size_t isam(0); isam < nCellsSamp.size(); ++isam ) {
1045 if ( isam == (size_t)CaloCell_ID::EME2 && std::get<1>(nCellsSamp.at(isam)) > 0 ) {
1047 }
1048 }
1049 nCellsSamp.clear();
1050 }
1051 }
1052
1053 return StatusCode::SUCCESS;
1054}
Scalar deltaPhi(const MatrixBase< Derived > &vec) const
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