11 #include "CLHEP/Geometry/Transform3D.h"
21 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
22 if(v_id[
i].
first.rawId()==
id.rawId()){
35 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
42 return std::pair<DetId, float>(id,
max);
49 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
56 return std::pair<DetId, float>(id,
max);
62 return getMaximum( cluster.hitsAndFractions(), recHits );
67 return getMaximum( cluster.hitsAndFractions(), recHits, flagsexcl, severitiesexcl, sevLv );
73 if (
id ==
DetId(0) ) {
77 if ( it != recHits->
end() ) {
78 return (*it).energy();
90 if (
id ==
DetId(0) ) {
94 if ( it != recHits->
end() ) {
96 uint32_t rhFlag = (*it).recoFlag();
97 std::vector<int>::const_iterator vit =
std::find( flagsexcl.begin(), flagsexcl.end(), rhFlag );
100 if ( vit != flagsexcl.end() )
return 0;
102 int severityFlag = sevLv->
severityLevel( it->id(), *recHits);
103 std::vector<int>::const_iterator sit =
std::find(severitiesexcl.begin(), severitiesexcl.end(), severityFlag);
106 if (sit!= severitiesexcl.end())
109 return (*it).energy();
138 const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
139 for (
int i = ixMin;
i <= ixMax; ++
i ) {
140 for (
int j = iyMin;
j <= iyMax; ++
j ) {
142 cursor.offsetBy(
i,
j );
156 float EcalClusterTools::matrixEnergy(
const reco::BasicCluster &cluster,
const EcalRecHitCollection *recHits,
const CaloTopology* topology,
DetId id,
int ixMin,
int ixMax,
int iyMin,
int iyMax,
const std::vector<int>& flagsexcl,
const std::vector<int>& severitiesexcl,
const EcalSeverityLevelAlgo *sevLv )
160 const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
162 for (
int i = ixMin;
i <= ixMax; ++
i ) {
163 for (
int j = iyMin;
j <= iyMax; ++
j ) {
165 cursor.offsetBy(
i,
j );
167 energy +=
recHitEnergy( *cursor, recHits, flagsexcl, severitiesexcl, sevLv )*
frac;
178 std::vector<DetId>
v;
179 for (
int i = ixMin;
i <= ixMax; ++
i ) {
180 for (
int j = iyMin;
j <= iyMax; ++
j ) {
182 cursor.offsetBy(
i,
j );
183 if ( *cursor !=
DetId(0) ) v.push_back( *cursor );
194 std::list<float> energies;
195 float max_E =
matrixEnergy( cluster, recHits, topology,
id, -1, 0, -1, 0 );
206 float max_E =
matrixEnergy( cluster, recHits, topology,
id, -1, 0, -1, 0,flagsexcl, severitiesexcl, sevLv );
207 max_E =
std::max( max_E,
matrixEnergy( cluster, recHits, topology,
id, -1, 0, 0, 1,flagsexcl, severitiesexcl, sevLv ) );
208 max_E =
std::max( max_E,
matrixEnergy( cluster, recHits, topology,
id, 0, 1, 0, 1,flagsexcl, severitiesexcl, sevLv ) );
209 max_E =
std::max( max_E,
matrixEnergy( cluster, recHits, topology,
id, 0, 1, -1, 0,flagsexcl, severitiesexcl, sevLv ) );
217 float max_E =
matrixEnergy( cluster, recHits, topology,
id, -1, 1, -1, 0 );
227 std::list<float> energies;
228 float max_E =
matrixEnergy( cluster, recHits, topology,
id, -1, 1, -1, 0,flagsexcl, severitiesexcl, sevLv );
229 max_E =
std::max( max_E,
matrixEnergy( cluster, recHits, topology,
id, 0, 1, -1, 1,flagsexcl, severitiesexcl, sevLv ) );
230 max_E =
std::max( max_E,
matrixEnergy( cluster, recHits, topology,
id, -1, 1, 0, 1,flagsexcl, severitiesexcl, sevLv ) );
231 max_E =
std::max( max_E,
matrixEnergy( cluster, recHits, topology,
id, -1, 0, -1, 1,flagsexcl, severitiesexcl, sevLv ) );
238 return matrixEnergy( cluster, recHits, topology,
id, -1, 1, -1, 1 );
245 return matrixEnergy( cluster, recHits, topology,
id, -1, 1, -1, 1,flagsexcl, severitiesexcl, sevLv );
252 float max_E =
matrixEnergy( cluster, recHits, topology,
id, -1, 2, -2, 1 );
262 std::list<float> energies;
263 float max_E =
matrixEnergy( cluster, recHits, topology,
id, -1, 2, -2, 1,flagsexcl, severitiesexcl, sevLv );
264 max_E =
std::max( max_E,
matrixEnergy( cluster, recHits, topology,
id, -2, 1, -2, 1,flagsexcl, severitiesexcl, sevLv ) );
265 max_E =
std::max( max_E,
matrixEnergy( cluster, recHits, topology,
id, -2, 1, -1, 2,flagsexcl, severitiesexcl, sevLv ) );
266 max_E =
std::max( max_E,
matrixEnergy( cluster, recHits, topology,
id, -1, 2, -1, 2,flagsexcl, severitiesexcl, sevLv ) );
275 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, -2, 2 );
280 DetId id =
getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).
first;
281 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, -2, 2,flagsexcl, severitiesexcl, sevLv );
291 return getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).
second;
297 std::vector<float> energies;
298 const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
299 energies.reserve( v_id.size() );
300 if ( v_id.size() < 2 )
return 0;
301 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
304 std::partial_sort( energies.begin(), energies.begin()+2, energies.end(), std::greater<float>() );
312 std::vector<float> energies;
313 const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
314 energies.reserve( v_id.size() );
315 if ( v_id.size() < 2 )
return 0;
316 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
319 std::partial_sort( energies.begin(), energies.begin()+2, energies.end(), std::greater<float>() );
330 return matrixEnergy( cluster, recHits, topology,
id, 1, 2, -2, 2 );
336 DetId id =
getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).
first;
337 return matrixEnergy( cluster, recHits, topology,
id, 1, 2, -2, 2,flagsexcl, severitiesexcl, sevLv );
344 return matrixEnergy( cluster, recHits, topology,
id, -2, -1, -2, 2 );
349 DetId id =
getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).
first;
350 return matrixEnergy( cluster, recHits, topology,
id, -2, -1, -2, 2,flagsexcl, severitiesexcl, sevLv );
358 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, 1, 2 );
363 DetId id =
getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).
first;
364 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, 1, 2,flagsexcl, severitiesexcl, sevLv );
371 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, -2, -1 );
376 DetId id =
getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).
first;
377 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, -2, -1,flagsexcl, severitiesexcl, sevLv );
387 float left =
matrixEnergy( cluster, recHits, topology,
id, -1, -1, -2, 2 );
389 float right =
matrixEnergy( cluster, recHits, topology,
id, 1, 1, -2, 2 );
391 float centre =
matrixEnergy( cluster, recHits, topology,
id, 0, 0, -2, 2 );
394 return left > right ? left+centre : right+centre;
399 DetId id =
getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).
first;
402 float left =
matrixEnergy( cluster, recHits, topology,
id, -1, -1, -2, 2,flagsexcl, severitiesexcl, sevLv );
404 float right =
matrixEnergy( cluster, recHits, topology,
id, 1, 1, -2, 2,flagsexcl, severitiesexcl, sevLv );
406 float centre =
matrixEnergy( cluster, recHits, topology,
id, 0, 0, -2, 2,flagsexcl, severitiesexcl, sevLv );
409 return left > right ? left+centre : right+centre;
416 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, -2, 2 );
422 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, -2, 2 ,
423 flagsexcl, severitiesexcl, sevLv);
431 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, 0, 0 );
437 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, 0, 0,flagsexcl, severitiesexcl, sevLv );
444 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, -1, 1 );
450 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, -1, 1,flagsexcl, severitiesexcl, sevLv );
457 return matrixEnergy( cluster, recHits, topology,
id, -1, 1, 0, 0 );
463 return matrixEnergy( cluster, recHits, topology,
id, -1, 1, 0, 0,flagsexcl, severitiesexcl, sevLv );
470 return matrixEnergy( cluster, recHits, topology,
id, -1, -1, 0, 0 );
476 return matrixEnergy( cluster, recHits, topology,
id, -1, -1, 0, 0,flagsexcl, severitiesexcl, sevLv );
483 return matrixEnergy( cluster, recHits, topology,
id, 1, 1, 0, 0 );
489 return matrixEnergy( cluster, recHits, topology,
id, 1, 1, 0, 0,flagsexcl, severitiesexcl, sevLv );
496 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, 1, 1 );
502 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, 1, 1,flagsexcl, severitiesexcl, sevLv );
510 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, -1, -1 );
516 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, -1, -1,flagsexcl, severitiesexcl, sevLv );
522 float clusterEnergy = cluster.energy();
523 const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
525 edm::LogWarning(
"EcalClusterTools::energyBasketFractionEta") <<
"Trying to get basket fraction for endcap basic-clusters. Basket fractions can be obtained ONLY for barrel basic-clusters. Returning empty vector.";
526 return basketFraction;
528 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
531 std::sort( basketFraction.rbegin(), basketFraction.rend() );
532 return basketFraction;
539 float clusterEnergy = cluster.energy();
540 const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
542 edm::LogWarning(
"EcalClusterTools::energyBasketFractionEta") <<
"Trying to get basket fraction for endcap basic-clusters. Basket fractions can be obtained ONLY for barrel basic-clusters. Returning empty vector.";
543 return basketFraction;
545 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
548 std::sort( basketFraction.rbegin(), basketFraction.rend() );
549 return basketFraction;
555 float clusterEnergy = cluster.energy();
556 const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
558 edm::LogWarning(
"EcalClusterTools::energyBasketFractionPhi") <<
"Trying to get basket fraction for endcap basic-clusters. Basket fractions can be obtained ONLY for barrel basic-clusters. Returning empty vector.";
559 return basketFraction;
561 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
564 std::sort( basketFraction.rbegin(), basketFraction.rend() );
565 return basketFraction;
572 float clusterEnergy = cluster.energy();
573 const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
575 edm::LogWarning(
"EcalClusterTools::energyBasketFractionPhi") <<
"Trying to get basket fraction for endcap basic-clusters. Basket fractions can be obtained ONLY for barrel basic-clusters. Returning empty vector.";
576 return basketFraction;
578 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
581 std::sort( basketFraction.rbegin(), basketFraction.rend() );
582 return basketFraction;
588 std::vector<EcalClusterTools::EcalClusterEnergyDeposition> energyDistribution;
591 CLHEP::Hep3Vector clVect(cluster.position().x(), cluster.position().y(), cluster.position().z());
592 CLHEP::Hep3Vector clDir(clVect);
593 clDir*=1.0/clDir.mag();
595 CLHEP::Hep3Vector theta_axis(clDir.y(),-clDir.x(),0.0);
596 theta_axis *= 1.0/theta_axis.mag();
597 CLHEP::Hep3Vector phi_axis = theta_axis.cross(clDir);
599 const std::vector< std::pair<DetId, float> >& clusterDetIds = cluster.hitsAndFractions();
603 std::vector< std::pair<DetId, float> >::const_iterator posCurrent;
605 for(posCurrent=clusterDetIds.begin(); posCurrent!=clusterDetIds.end(); ++posCurrent) {
609 if(( (*posCurrent).first !=
DetId(0)) && (recHits->
find( (*posCurrent).first ) != recHits->
end())) {
617 LogDebug(
"ClusterShapeAlgo") <<
"Crystal has insufficient energy: E = "
623 DetId id_ = (*posCurrent).first;
626 CLHEP::Hep3Vector gblPos (cellPos.
x(),cellPos.
y(),cellPos.
z());
628 CLHEP::Hep3Vector
diff = gblPos - clVect;
632 CLHEP::Hep3Vector DigiVect = diff - diff.dot(clDir)*clDir;
633 clEdep.
r = DigiVect.mag();
635 <<
"\tdiff = " << diff.mag()
636 <<
"\tr = " << clEdep.
r;
637 clEdep.
phi = DigiVect.angle(theta_axis);
638 if(DigiVect.dot(phi_axis)<0) clEdep.
phi = 2 *
M_PI - clEdep.
phi;
639 energyDistribution.push_back(clEdep);
642 return energyDistribution;
649 std::vector<EcalClusterTools::EcalClusterEnergyDeposition> energyDistribution =
getEnergyDepTopology( cluster, recHits, geometry, logW, w0 );
651 std::vector<float>
lat;
652 double r, redmoment=0;
653 double phiRedmoment = 0 ;
654 double etaRedmoment = 0 ;
656 int clusterSize=energyDistribution.size();
657 float etaLat_, phiLat_, lat_;
668 if (energyDistribution[1].deposited_energy >
669 energyDistribution[0].deposited_energy)
671 tmp=n2; n2=n1; n1=
tmp;
673 for (
int i=2;
i<clusterSize;
i++) {
675 if (energyDistribution[
i].deposited_energy >
676 energyDistribution[n1].deposited_energy)
679 n2 = n1; n1 =
i; n=
tmp;
681 if (energyDistribution[
i].deposited_energy >
682 energyDistribution[n2].deposited_energy)
688 r = energyDistribution[
n].r;
689 redmoment += r*r* energyDistribution[
n].deposited_energy;
690 double rphi = r *
cos (energyDistribution[n].
phi) ;
691 phiRedmoment += rphi * rphi * energyDistribution[
n].deposited_energy;
692 double reta = r *
sin (energyDistribution[n].phi) ;
693 etaRedmoment += reta * reta * energyDistribution[
n].deposited_energy;
695 double e1 = energyDistribution[n1].deposited_energy;
696 double e2 = energyDistribution[n2].deposited_energy;
698 lat_ = redmoment/(redmoment+2.19*2.19*(e1+
e2));
699 phiLat_ = phiRedmoment/(phiRedmoment+2.19*2.19*(e1+
e2));
700 etaLat_ = etaRedmoment/(etaRedmoment+2.19*2.19*(e1+
e2));
702 lat.push_back(etaLat_);
703 lat.push_back(phiLat_);
714 const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
716 for(
const std::pair<DetId,float>& hitAndFrac : hsAndFs ) {
717 for( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
718 if( hitAndFrac.first != *it )
continue;
724 return meanPosition /
e5x5( cluster, recHits, topology );
733 const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
735 for(
const std::pair<DetId,float>& hitAndFrac : hsAndFs ) {
736 for( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
737 if( hitAndFrac.first != *it )
continue;
740 meanPosition = meanPosition +
recHitEnergy( *it, recHits,flagsexcl, severitiesexcl, sevLv ) *
position * hitAndFrac.second;
743 return meanPosition /
e5x5( cluster, recHits, topology,flagsexcl, severitiesexcl, sevLv );
761 const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
762 std::vector<DetId> v_id =
matrixDetId( topology,seedId, -2, 2, -2, 2 );
763 for(
const std::pair<DetId,float>& hAndF : hsAndFs ) {
764 for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
765 if( hAndF.first != *it )
continue;
767 if(energy<0.)
continue;
775 return std::pair<float,float>(meanDEta,meanDPhi);
785 const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
786 std::vector<DetId> v_id =
matrixDetId( topology,seedId, -2, 2, -2, 2 );
787 for(
const std::pair<DetId,float>& hAndF : hsAndFs ) {
788 for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
789 if( hAndF.first != *it )
continue;
790 float energy =
recHitEnergy(*it,recHits,flagsexcl, severitiesexcl, sevLv) * hAndF.second;
791 if(energy<0.)
continue;
799 return std::pair<float,float>(meanDEta,meanDPhi);
812 std::pair<float,float> meanXY(0.,0.);
817 const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
818 std::vector<DetId> v_id =
matrixDetId( topology,seedId, -2, 2, -2, 2 );
819 for(
const std::pair<DetId,float>& hAndF : hsAndFs ) {
820 for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
821 if( hAndF.first != *it )
continue;
823 if(energy<0.)
continue;
838 std::pair<float,float> meanXY(0.,0.);
843 const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
844 std::vector<DetId> v_id =
matrixDetId( topology,seedId, -2, 2, -2, 2 );
845 for(
const std::pair<DetId,float>& hAndF : hsAndFs ) {
846 for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
847 if( hAndF.first != *it )
continue;
848 float energy =
recHitEnergy(*it,recHits,flagsexcl, severitiesexcl, sevLv) * hAndF.second;
849 if(energy<0.)
continue;
863 float e_5x5 =
e5x5( cluster, recHits, topology );
864 float covEtaEta, covEtaPhi, covPhiPhi;
867 const std::vector< std::pair<DetId, float>>& v_id =cluster.hitsAndFractions();
871 double numeratorEtaEta = 0;
872 double numeratorEtaPhi = 0;
873 double numeratorPhiPhi = 0;
878 for (
int i = -2;
i <= 2; ++
i ) {
879 for (
int j = -2;
j <= 2; ++
j ) {
881 cursor.offsetBy(
i,
j );
885 if ( energy <= 0 )
continue;
889 double dPhi = position.
phi() - meanPosition.
phi();
893 double dEta = position.
eta() - meanPosition.eta();
898 numeratorEtaEta += w * dEta * dEta;
899 numeratorEtaPhi += w * dEta *
dPhi;
900 numeratorPhiPhi += w * dPhi *
dPhi;
904 if (denominator != 0.0) {
921 std::vector<float>
v;
922 v.push_back( covEtaEta );
923 v.push_back( covEtaPhi );
924 v.push_back( covPhiPhi );
932 float e_5x5 =
e5x5( cluster, recHits, topology,flagsexcl, severitiesexcl, sevLv );
933 float covEtaEta, covEtaPhi, covPhiPhi;
936 const std::vector<std::pair<DetId, float>>& v_id= cluster.hitsAndFractions();
940 double numeratorEtaEta = 0;
941 double numeratorEtaPhi = 0;
942 double numeratorPhiPhi = 0;
948 for (
int i = -2;
i <= 2; ++
i ) {
949 for (
int j = -2;
j <= 2; ++
j ) {
951 cursor.offsetBy(
i,
j );
955 if ( energy <= 0 )
continue;
959 double dPhi = position.
phi() - meanPosition.
phi();
963 double dEta = position.
eta() - meanPosition.eta();
968 numeratorEtaEta += w * dEta * dEta;
969 numeratorEtaPhi += w * dEta *
dPhi;
970 numeratorPhiPhi += w * dPhi *
dPhi;
974 if (denominator != 0.0) {
991 std::vector<float>
v;
992 v.push_back( covEtaEta );
993 v.push_back( covEtaPhi );
994 v.push_back( covPhiPhi );
1009 float e_5x5 =
e5x5( cluster, recHits, topology );
1010 float covEtaEta, covEtaPhi, covPhiPhi;
1014 const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
1019 double numeratorEtaEta = 0;
1020 double numeratorEtaPhi = 0;
1021 double numeratorPhiPhi = 0;
1026 const double barrelCrysSize = 0.01745;
1027 const double endcapCrysSize = 0.0447;
1032 const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1034 CaloNavigator<DetId> cursor = CaloNavigator<DetId>( seedId, topology->
getSubdetectorTopology( seedId ) );
1036 for (
int eastNr = -2; eastNr <= 2; ++eastNr ) {
1037 for (
int northNr = -2; northNr <= 2; ++northNr ) {
1039 cursor.offsetBy( eastNr, northNr);
1042 if ( energy <= 0 )
continue;
1047 if(isBarrel) dPhi =
getNrCrysDiffInPhi(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.second;
1048 else dPhi =
getDPhiEndcap(*cursor,mean5x5XYPos.first,mean5x5XYPos.second);
1054 numeratorEtaEta += w * dEta * dEta;
1055 numeratorEtaPhi += w * dEta *
dPhi;
1056 numeratorPhiPhi += w * dPhi *
dPhi;
1062 if (denominator != 0.0) {
1063 covEtaEta = crysSize*crysSize* numeratorEtaEta /
denominator;
1064 covEtaPhi = crysSize*crysSize* numeratorEtaPhi /
denominator;
1065 covPhiPhi = crysSize*crysSize* numeratorPhiPhi /
denominator;
1080 std::vector<float>
v;
1081 v.push_back( covEtaEta );
1082 v.push_back( covEtaPhi );
1083 v.push_back( covPhiPhi );
1092 float e_5x5 =
e5x5( cluster, recHits, topology,flagsexcl, severitiesexcl, sevLv );
1093 float covEtaEta, covEtaPhi, covPhiPhi;
1097 const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
1099 std::pair<float,float> mean5x5XYPos =
mean5x5PositionInXY(cluster,recHits,topology,flagsexcl, severitiesexcl, sevLv);
1102 double numeratorEtaEta = 0;
1103 double numeratorEtaPhi = 0;
1104 double numeratorPhiPhi = 0;
1109 const double barrelCrysSize = 0.01745;
1110 const double endcapCrysSize = 0.0447;
1115 const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1117 CaloNavigator<DetId> cursor = CaloNavigator<DetId>( seedId, topology->
getSubdetectorTopology( seedId ) );
1119 for (
int eastNr = -2; eastNr <= 2; ++eastNr ) {
1120 for (
int northNr = -2; northNr <= 2; ++northNr ) {
1122 cursor.offsetBy( eastNr, northNr);
1125 if ( energy <= 0 )
continue;
1130 if(isBarrel) dPhi =
getNrCrysDiffInPhi(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.second;
1131 else dPhi =
getDPhiEndcap(*cursor,mean5x5XYPos.first,mean5x5XYPos.second);
1137 numeratorEtaEta += w * dEta * dEta;
1138 numeratorEtaPhi += w * dEta *
dPhi;
1139 numeratorPhiPhi += w * dPhi *
dPhi;
1145 if (denominator != 0.0) {
1146 covEtaEta = crysSize*crysSize* numeratorEtaEta /
denominator;
1147 covEtaPhi = crysSize*crysSize* numeratorEtaPhi /
denominator;
1148 covPhiPhi = crysSize*crysSize* numeratorPhiPhi /
denominator;
1163 std::vector<float>
v;
1164 v.push_back( covEtaEta );
1165 v.push_back( covEtaPhi );
1166 v.push_back( covPhiPhi );
1188 if ((m>n) || ((n-m)%2 != 0) || (n<0) || (m<0))
return -1;
1192 if ((n>20) || (R0<=2.19))
return -1;
1200 double r,ph,
e,Re=0,Im=0;
1201 double TotalEnergy = cluster.energy();
1202 int index = (n/2)*(n/2)+(n/2)+m;
1203 std::vector<EcalClusterEnergyDeposition> energyDistribution =
getEnergyDepTopology( cluster, recHits, geometry, logW, w0 );
1204 int clusterSize = energyDistribution.size();
1205 if(clusterSize < 3)
return 0.0;
1207 for (
int i=0;
i<clusterSize;
i++)
1209 r = energyDistribution[
i].r / R0;
1211 std::vector<double> pol;
1212 pol.push_back(
f00(r) );
1213 pol.push_back(
f11(r) );
1214 pol.push_back(
f20(r) );
1215 pol.push_back(
f22(r) );
1216 pol.push_back(
f31(r) );
1217 pol.push_back(
f33(r) );
1218 pol.push_back(
f40(r) );
1219 pol.push_back(
f42(r) );
1220 pol.push_back(
f44(r) );
1221 pol.push_back(
f51(r) );
1222 pol.push_back(
f53(r) );
1223 pol.push_back(
f55(r) );
1224 ph = (energyDistribution[
i]).
phi;
1225 e = energyDistribution[
i].deposited_energy;
1226 Re = Re + e/TotalEnergy * pol[
index] *
cos( (
double) m * ph);
1227 Im = Im - e/TotalEnergy * pol[
index] *
sin( (
double) m * ph);
1230 return sqrt(Re*Re+Im*Im);
1237 double r, ph,
e, Re=0, Im=0, f_nm;
1238 double TotalEnergy = cluster.energy();
1239 std::vector<EcalClusterEnergyDeposition> energyDistribution =
getEnergyDepTopology( cluster, recHits, geometry, logW, w0 );
1240 int clusterSize=energyDistribution.size();
1241 if(clusterSize<3)
return 0.0;
1243 for (
int i = 0;
i < clusterSize; ++
i)
1245 r = energyDistribution[
i].r / R0;
1247 ph = energyDistribution[
i].phi;
1248 e = energyDistribution[
i].deposited_energy;
1250 for (
int s=0;
s<=(n-
m)/2;
s++) {
1257 Re = Re + e/TotalEnergy * f_nm *
cos( (
double) m*ph);
1258 Im = Im - e/TotalEnergy * f_nm *
sin( (
double) m*ph);
1261 return sqrt(Re*Re+Im*Im);
1278 return std::sqrt(iXNorm*iXNorm+iYNorm*iYNorm);
1305 int iXNorm = eeId.
ix()-50;
1306 if(iXNorm<=0) iXNorm--;
1317 int iYNorm = eeId.
iy()-50;
1318 if(iYNorm<=0) iYNorm--;
1327 float crysIEta =
getIEta(crysId);
1328 float orginIEta =
getIEta(orginId);
1331 float nrCrysDiff = crysIEta-orginIEta;
1336 if(crysIEta*orginIEta<0){
1337 if(crysIEta>0) nrCrysDiff--;
1347 float crysIPhi =
getIPhi(crysId);
1348 float orginIPhi =
getIPhi(orginId);
1351 float nrCrysDiff = crysIPhi-orginIPhi;
1354 if (nrCrysDiff > + 180) { nrCrysDiff = nrCrysDiff - 360; }
1355 if (nrCrysDiff < - 180) { nrCrysDiff = nrCrysDiff + 360; }
1366 float hitLocalR2 = (iXNorm-meanX)*(iXNorm-meanX)+(iYNorm-meanY)*(iYNorm-meanY);
1367 float hitR2 = iXNorm*iXNorm+iYNorm*iYNorm;
1368 float meanR2 = meanX*meanX+meanY*meanY;
1369 float hitR =
sqrt(hitR2);
1370 float meanR =
sqrt(meanR2);
1372 float tmp = (hitR2+meanR2-hitLocalR2)/(2*hitR*meanR);
1373 if (tmp<-1) tmp =-1;
1375 float phi = acos(tmp);
1385 float e_5x5 =
e5x5(bcluster, recHits, topology);
1386 float covEtaEta, covEtaPhi, covPhiPhi;
1389 const std::vector<std::pair<DetId, float> >& v_id = cluster.
hitsAndFractions();
1393 double numeratorEtaEta = 0;
1394 double numeratorEtaPhi = 0;
1395 double numeratorPhiPhi = 0;
1398 const double barrelCrysSize = 0.01745;
1399 const double endcapCrysSize = 0.0447;
1404 const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1406 for (
size_t i = 0;
i < v_id.size(); ++
i) {
1411 if (energy <= 0)
continue;
1415 if(isBarrel) dPhi =
getNrCrysDiffInPhi(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.second;
1416 else dPhi =
getDPhiEndcap(*cursor,mean5x5XYPos.first,mean5x5XYPos.second);
1424 numeratorEtaEta += w * dEta * dEta;
1425 numeratorEtaPhi += w * dEta *
dPhi;
1426 numeratorPhiPhi += w * dPhi *
dPhi;
1430 if (denominator != 0.0) {
1431 covEtaEta = crysSize*crysSize* numeratorEtaEta /
denominator;
1432 covEtaPhi = crysSize*crysSize* numeratorEtaPhi /
denominator;
1433 covPhiPhi = crysSize*crysSize* numeratorPhiPhi /
denominator;
1448 std::vector<float>
v;
1449 v.push_back( covEtaEta );
1450 v.push_back( covEtaPhi );
1451 v.push_back( covPhiPhi );
1463 float e_5x5 =
e5x5(bcluster, recHits, topology);
1464 float covEtaEta, covEtaPhi, covPhiPhi;
1467 const std::vector<std::pair<DetId, float> >& v_id = cluster.
hitsAndFractions();
1469 std::pair<float,float> mean5x5XYPos =
mean5x5PositionInXY(cluster,recHits,topology,flagsexcl, severitiesexcl, sevLv);
1471 double numeratorEtaEta = 0;
1472 double numeratorEtaPhi = 0;
1473 double numeratorPhiPhi = 0;
1476 const double barrelCrysSize = 0.01745;
1477 const double endcapCrysSize = 0.0447;
1482 const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1484 for (
size_t i = 0;
i < v_id.size(); ++
i) {
1489 if (energy <= 0)
continue;
1493 if(isBarrel) dPhi =
getNrCrysDiffInPhi(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.second;
1494 else dPhi =
getDPhiEndcap(*cursor,mean5x5XYPos.first,mean5x5XYPos.second);
1502 numeratorEtaEta += w * dEta * dEta;
1503 numeratorEtaPhi += w * dEta *
dPhi;
1504 numeratorPhiPhi += w * dPhi *
dPhi;
1508 if (denominator != 0.0) {
1509 covEtaEta = crysSize*crysSize* numeratorEtaEta /
denominator;
1510 covEtaPhi = crysSize*crysSize* numeratorEtaPhi /
denominator;
1511 covPhiPhi = crysSize*crysSize* numeratorPhiPhi /
denominator;
1526 std::vector<float>
v;
1527 v.push_back( covEtaEta );
1528 v.push_back( covEtaPhi );
1529 v.push_back( covPhiPhi );
1546 std::vector<std::pair<const EcalRecHit*, float> > RH_ptrs_fracs;
1548 const std::vector< std::pair<DetId, float> >& myHitsPair = basicCluster.hitsAndFractions();
1550 for(
unsigned int i=0;
i<myHitsPair.size();
i++){
1553 RH_ptrs_fracs.push_back( std::make_pair(&(*myRH) , myHitsPair[
i].
second) );
1564 returnMoments.
sMaj = -1.;
1565 returnMoments.
sMin = -1.;
1566 returnMoments.
alpha = 0.;
1573 return returnMoments;
1580 double mid_eta(0),mid_phi(0),mid_x(0),mid_y(0);
1584 double max_phi=-10.;
1585 double min_phi=100.;
1588 std::vector<double> etaDetId;
1589 std::vector<double> phiDetId;
1590 std::vector<double> xDetId;
1591 std::vector<double> yDetId;
1592 std::vector<double> wiDetId;
1594 unsigned int nCry=0;
1599 for(std::vector<std::pair<const EcalRecHit*, float> >::const_iterator rhf_ptr = RH_ptrs_fracs.begin(); rhf_ptr != RH_ptrs_fracs.end(); rhf_ptr++){
1604 double temp_eta(0),temp_phi(0),temp_x(0),temp_y(0);
1617 double temp_ene=rh_ptr->
energy() * rhf_ptr->second;
1619 double temp_wi=((useLogWeights) ?
1624 if(temp_phi>max_phi) max_phi=temp_phi;
1625 if(temp_phi<min_phi) min_phi=temp_phi;
1626 etaDetId.push_back(temp_eta);
1627 phiDetId.push_back(temp_phi);
1628 xDetId.push_back(temp_x);
1629 yDetId.push_back(temp_y);
1630 wiDetId.push_back(temp_wi);
1631 denominator+=temp_wi;
1637 if(max_phi==359.5 && min_phi==0.5){
1638 for(
unsigned int i=0;
i<nCry;
i++){
1639 if(phiDetId[
i] - 179. > 0.) phiDetId[
i]-=360.;
1640 mid_phi+=phiDetId[
i]*wiDetId[
i];
1641 mid_eta+=etaDetId[
i]*wiDetId[
i];
1644 for(
unsigned int i=0;
i<nCry;
i++){
1645 mid_phi+=phiDetId[
i]*wiDetId[
i];
1646 mid_eta+=etaDetId[
i]*wiDetId[
i];
1650 for(
unsigned int i=0;
i<nCry;
i++){
1651 mid_eta+=etaDetId[
i]*wiDetId[
i];
1652 mid_x+=xDetId[
i]*wiDetId[
i];
1653 mid_y+=yDetId[
i]*wiDetId[
i];
1669 double deta(0),dphi(0);
1671 for(
unsigned int i=0;
i<nCry;
i++) {
1673 deta = etaDetId[
i]-mid_eta;
1674 dphi = phiDetId[
i]-mid_phi;
1676 deta = etaDetId[
i]-mid_eta;
1677 float hitLocalR2 = (xDetId[
i]-mid_x)*(xDetId[
i]-mid_x)+(yDetId[
i]-mid_y)*(yDetId[
i]-mid_y);
1678 float hitR2 = xDetId[
i]*xDetId[
i]+yDetId[
i]*yDetId[
i];
1679 float meanR2 = mid_x*mid_x+mid_y*mid_y;
1680 float hitR =
sqrt(hitR2);
1681 float meanR =
sqrt(meanR2);
1682 float phi = acos((hitR2+meanR2-hitLocalR2)/(2*hitR*meanR));
1686 See += (wiDetId[
i]* deta * deta) / denominator;
1687 Spp += phiCorrectionFactor*(wiDetId[
i]* dphi * dphi) / denominator;
1688 Sep +=
sqrt(phiCorrectionFactor)*(wiDetId[
i]*deta*dphi) / denominator;
1694 returnMoments.
sMaj = ((See + Spp) +
sqrt((See - Spp)*(See - Spp) + 4.*Sep*Sep)) / 2.;
1695 returnMoments.
sMin = ((See + Spp) -
sqrt((See - Spp)*(See - Spp) + 4.*Sep*Sep)) / 2.;
1697 returnMoments.
alpha = atan( (See - Spp +
sqrt( (Spp - See)*(Spp - See) + 4.*Sep*Sep )) / (2.*Sep));
1699 return returnMoments;
1714 std::vector<std::pair<const EcalRecHit*, float> > RH_ptrs_fracs;
1715 const std::vector< std::pair<DetId, float> >& myHitsPair = superCluster.
hitsAndFractions();
1716 for(
unsigned int i=0;
i< myHitsPair.size(); ++
i){
1719 if( myRH != recHits.
end() && myRH->energy()*myHitsPair[
i].second > energyThreshold){
1720 RH_ptrs_fracs.push_back( std::make_pair(&(*myRH) , myHitsPair[
i].
second) );
1734 std::vector<std::pair<const EcalRecHit*, float> > RH_ptrs_fracs;
1735 const std::vector< std::pair<DetId, float> >& myHitsPair = superCluster.
hitsAndFractions();
1736 for(
unsigned int i=0;
i<myHitsPair.size(); ++
i){
1739 if(myRH != recHits.
end() && myRH->energy()*myHitsPair[
i].second > energyRHThresh)
1740 RH_ptrs_fracs.push_back( std::make_pair(&(*myRH) , myHitsPair[
i].second) );
1747 EBDetId EBdetIdi( rh->detid() );
1748 float the_fraction = 0;
1750 bool inEtaWindow = (
abs(
deltaIEta(seedPosition[0],EBdetIdi.ieta()) ) <= ieta_delta );
1751 bool inPhiWindow = (
abs(
deltaIPhi(seedPosition[1],EBdetIdi.iphi()) ) <= iphi_delta );
1752 bool passEThresh = ( rh->energy() > energyRHThresh );
1753 bool alreadyCounted =
false;
1756 bool is_SCrh_inside_recHits =
false;
1757 for(
unsigned int i=0;
i<myHitsPair.size();
i++){
1759 if(SCrh != recHits.
end()){
1760 the_fraction = myHitsPair[
i].second;
1761 is_SCrh_inside_recHits =
true;
1762 if( rh->detid() == SCrh->detid() ) alreadyCounted =
true;
1766 if( is_SCrh_inside_recHits && !alreadyCounted && passEThresh && inEtaWindow && inPhiWindow){
1767 RH_ptrs_fracs.push_back( std::make_pair(&(*rh),the_fraction) );
1780 std::vector<float> shapes;
1783 if(RH_ptrs_fracs.size()<2){
1784 shapes.push_back( -3 );
1785 shapes.push_back( -3 );
1791 int tempInt = seedPosition[0];
1792 if(tempInt <0) tempInt++;
1796 float centerIEta = 0.;
1797 float centerIPhi = 0.;
1800 for(std::vector<std::pair<const EcalRecHit*,float> >::const_iterator rhf_ptr = RH_ptrs_fracs.begin(); rhf_ptr != RH_ptrs_fracs.end(); rhf_ptr++){
1804 if(fabs(energyTotal) < 0.0001){
1806 shapes.push_back( -2 );
1807 shapes.push_back( -2 );
1810 float rh_energy = rh_ptr->
energy() * rhf_ptr->second;
1812 if(fabs(weightedPositionMethod)<0.0001){
1813 weight = rh_energy/energyTotal;
1815 weight =
std::max(0.0, 4.2 +
log(rh_energy/energyTotal));
1818 centerIEta += weight*
deltaIEta(seedPosition[0],EBdetIdi.ieta());
1819 centerIPhi += weight*
deltaIPhi(seedPosition[1],EBdetIdi.iphi());
1821 if(fabs(denominator) < 0.0001){
1823 shapes.push_back( -2 );
1824 shapes.push_back( -2 );
1832 TMatrixDSym inertia(2);
1833 double inertia00 = 0.;
1834 double inertia01 = 0.;
1835 double inertia11 = 0.;
1837 for(std::vector<std::pair<const EcalRecHit*,float> >::const_iterator rhf_ptr = RH_ptrs_fracs.begin(); rhf_ptr != RH_ptrs_fracs.end(); rhf_ptr++){
1842 if(fabs(energyTotal) < 0.0001){
1844 shapes.push_back( -2 );
1845 shapes.push_back( -2 );
1848 float rh_energy = rh_ptr->
energy() * rhf_ptr->second;
1850 if(fabs(weightedPositionMethod) < 0.0001){
1851 weight = rh_energy/energyTotal;
1853 weight =
std::max(0.0, 4.2 +
log(rh_energy/energyTotal));
1856 float ieta_rh_to_center =
deltaIEta(seedPosition[0],EBdetIdi.ieta()) - centerIEta;
1857 float iphi_rh_to_center =
deltaIPhi(seedPosition[1],EBdetIdi.iphi()) - centerIPhi;
1859 inertia00 += weight*iphi_rh_to_center*iphi_rh_to_center;
1860 inertia01 -= weight*iphi_rh_to_center*ieta_rh_to_center;
1861 inertia11 += weight*ieta_rh_to_center*ieta_rh_to_center;
1865 inertia[0][0] = inertia00;
1866 inertia[0][1] = inertia01;
1867 inertia[1][0] = inertia01;
1868 inertia[1][1] = inertia11;
1872 TMatrixD eVectors(2,2);
1873 TVectorD eValues(2);
1875 eVectors=inertia.EigenVectors(eValues);
1882 TVectorD smallerAxis(2);
1883 smallerAxis[0]=eVectors[0][1];
1884 smallerAxis[1]=eVectors[1][1];
1887 Double_t
temp = fabs(smallerAxis[0]);
1888 if(fabs(eValues[0]) < 0.0001){
1890 shapes.push_back( -2 );
1891 shapes.push_back( -2 );
1895 float Roundness = eValues[1]/eValues[0];
1896 float Angle=acos(temp);
1898 if( -0.00001 < Roundness && Roundness < 0) Roundness = 0.;
1899 if( -0.00001 < Angle && Angle < 0 ) Angle = 0.;
1901 shapes.push_back( Roundness );
1902 shapes.push_back( Angle );
1911 int rel_iphi = rh_iphi - seed_iphi;
1913 if(rel_iphi > 180) rel_iphi = rel_iphi - 360;
1914 if(rel_iphi < -180) rel_iphi = rel_iphi + 360;
1923 if(seed_ieta < 0) seed_ieta++;
1924 if(rh_ieta < 0) rh_ieta++;
1925 int rel_ieta = rh_ieta - seed_ieta;
1931 std::vector<int> seedPosition;
1936 for(std::vector<std::pair<const EcalRecHit*,float> >::const_iterator rhf_ptr = RH_ptrs_fracs.begin(); rhf_ptr != RH_ptrs_fracs.end(); rhf_ptr++){
1940 float rh_energy = rh_ptr->
energy() * rhf_ptr->second;
1942 if(eSeedRH < rh_energy){
1943 eSeedRH = rh_energy;
1944 iEtaSeedRH = EBdetIdi.ieta();
1945 iPhiSeedRH = EBdetIdi.iphi();
1950 seedPosition.push_back(iEtaSeedRH);
1951 seedPosition.push_back(iPhiSeedRH);
1952 return seedPosition;
1958 for(
const auto& hAndF : RH_ptrs_fracs ) {
1959 sumE += hAndF.first->energy() * hAndF.second;
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
const DetId & detid() const
Sin< T >::type sin(const T &t)
Geom::Phi< T > phi() const
std::vector< EcalRecHit >::const_iterator const_iterator
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
static const int kTowersInPhi
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
static int position[TOTALCHAMBERS][3]
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
int iphi() const
get the crystal iphi
U second(std::pair< T, U > const &p)
int im() const
get the number of module inside the SM (1-4)
static const int kCrystalsInPhi
double dPhi(double phi1, double phi2)
const T & max(const T &a, const T &b)
Cos< T >::type cos(const T &t)
static const int kModulesPerSM
Abs< T >::type abs(const T &t)
int ieta() const
get the crystal ieta
int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
const_iterator end() const
static const int MAX_IPHI
XYZVectorD XYZVector
spatial vector with cartesian internal representation
const CaloSubdetectorTopology * getSubdetectorTopology(const DetId &id) const
access the subdetector Topology for the given subdetector directly
std::vector< std::vector< double > > tmp
ESHandle< TrackerGeometry > geometry
iterator find(key_type k)
const CaloClusterPtr & seed() const
seed BasicCluster
double energySum(const DataFrame &df, int fs, int ls)
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
EcalSeverityLevel::SeverityLevel severityLevel(const DetId &id, const EcalRecHitCollection &rhs) const
Evaluate status from id.
Power< A, B >::type pow(const A &a, const B &b)
const_iterator begin() const