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()){
34 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
41 return std::pair<DetId, float>(
id,
max);
48 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
55 return std::pair<DetId, float>(
id,
max);
61 return getMaximum( cluster.hitsAndFractions(), recHits );
66 return getMaximum( cluster.hitsAndFractions(), recHits, flagsexcl, severitiesexcl, sevLv );
72 if (
id ==
DetId(0) ) {
76 if ( it != recHits->
end() ) {
77 return (*it).energy();
89 if (
id ==
DetId(0) ) {
93 if ( it != recHits->
end() ) {
95 uint32_t rhFlag = (*it).recoFlag();
96 std::vector<int>::const_iterator vit =
std::find( flagsexcl.begin(), flagsexcl.end(), rhFlag );
99 if ( vit != flagsexcl.end() )
return 0;
101 int severityFlag = sevLv->
severityLevel( it->id(), *recHits);
102 std::vector<int>::const_iterator sit =
std::find(severitiesexcl.begin(), severitiesexcl.end(), severityFlag);
105 if (sit!= severitiesexcl.end())
108 return (*it).energy();
137 std::vector< std::pair<DetId, float> > v_id = cluster.hitsAndFractions();
138 for (
int i = ixMin;
i <= ixMax; ++
i ) {
139 for (
int j = iyMin;
j <= iyMax; ++
j ) {
155 float EcalClusterTools::matrixEnergy(
const reco::BasicCluster &cluster,
const EcalRecHitCollection *recHits,
const CaloTopology* topology,
DetId id,
int ixMin,
int ixMax,
int iyMin,
int iyMax, std::vector<int> flagsexcl, std::vector<int> severitiesexcl,
const EcalSeverityLevelAlgo *sevLv )
160 for (
int i = ixMin;
i <= ixMax; ++
i ) {
161 for (
int j = iyMin;
j <= iyMax; ++
j ) {
164 energy +=
recHitEnergy( *cursor, recHits, flagsexcl, severitiesexcl, sevLv );
175 std::vector<DetId>
v;
176 for (
int i = ixMin;
i <= ixMax; ++
i ) {
177 for (
int j = iyMin;
j <= iyMax; ++
j ) {
180 if ( *cursor !=
DetId(0) ) v.push_back( *cursor );
191 std::list<float> energies;
192 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -1, 0, -1, 0 ) );
193 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -1, 0, 0, 1 ) );
194 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, 0, 1, 0, 1 ) );
195 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, 0, 1, -1, 0 ) );
198 return *std::max_element(energies.begin(),energies.end());
206 std::list<float> energies;
207 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -1, 0, -1, 0,flagsexcl, severitiesexcl, sevLv ) );
208 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -1, 0, 0, 1,flagsexcl, severitiesexcl, sevLv ) );
209 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, 0, 1, 0, 1,flagsexcl, severitiesexcl, sevLv ) );
210 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, 0, 1, -1, 0,flagsexcl, severitiesexcl, sevLv ) );
213 return *std::max_element(energies.begin(),energies.end());
221 std::list<float> energies;
222 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -1, 1, -1, 0 ) );
223 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, 0, 1, -1, 1 ) );
224 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -1, 1, 0, 1 ) );
225 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -1, 0, -1, 1 ) );
226 return *std::max_element(energies.begin(),energies.end());
232 std::list<float> energies;
233 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -1, 1, -1, 0,flagsexcl, severitiesexcl, sevLv ) );
234 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, 0, 1, -1, 1,flagsexcl, severitiesexcl, sevLv ) );
235 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -1, 1, 0, 1,flagsexcl, severitiesexcl, sevLv ) );
236 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -1, 0, -1, 1,flagsexcl, severitiesexcl, sevLv ) );
237 return *std::max_element(energies.begin(),energies.end());
243 return matrixEnergy( cluster, recHits, topology,
id, -1, 1, -1, 1 );
250 return matrixEnergy( cluster, recHits, topology,
id, -1, 1, -1, 1,flagsexcl, severitiesexcl, sevLv );
257 std::list<float> energies;
258 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -1, 2, -2, 1 ) );
259 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -2, 1, -2, 1 ) );
260 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -2, 1, -1, 2 ) );
261 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -1, 2, -1, 2 ) );
262 return *std::max_element(energies.begin(),energies.end());
268 std::list<float> energies;
269 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -1, 2, -2, 1,flagsexcl, severitiesexcl, sevLv ) );
270 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -2, 1, -2, 1,flagsexcl, severitiesexcl, sevLv ) );
271 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -2, 1, -1, 2,flagsexcl, severitiesexcl, sevLv ) );
272 energies.push_back(
matrixEnergy( cluster, recHits, topology,
id, -1, 2, -1, 2,flagsexcl, severitiesexcl, sevLv ) );
273 return *std::max_element(energies.begin(),energies.end());
281 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, -2, 2 );
286 DetId id =
getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).
first;
287 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, -2, 2,flagsexcl, severitiesexcl, sevLv );
297 return getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).
second;
303 std::list<float> energies;
304 std::vector< std::pair<DetId, float> > v_id = cluster.hitsAndFractions();
305 if ( v_id.size() < 2 )
return 0;
306 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
310 return *--(--energies.end());
317 std::list<float> energies;
318 std::vector< std::pair<DetId, float> > v_id = cluster.hitsAndFractions();
319 if ( v_id.size() < 2 )
return 0;
320 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
324 return *--(--energies.end());
334 return matrixEnergy( cluster, recHits, topology,
id, 1, 2, -2, 2 );
340 DetId id =
getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).
first;
341 return matrixEnergy( cluster, recHits, topology,
id, 1, 2, -2, 2,flagsexcl, severitiesexcl, sevLv );
348 return matrixEnergy( cluster, recHits, topology,
id, -2, -1, -2, 2 );
353 DetId id =
getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).
first;
354 return matrixEnergy( cluster, recHits, topology,
id, -2, -1, -2, 2,flagsexcl, severitiesexcl, sevLv );
362 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, 1, 2 );
367 DetId id =
getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).
first;
368 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, 1, 2,flagsexcl, severitiesexcl, sevLv );
375 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, -2, -1 );
380 DetId id =
getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).
first;
381 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, -2, -1,flagsexcl, severitiesexcl, sevLv );
391 float left =
matrixEnergy( cluster, recHits, topology,
id, -1, -1, -2, 2 );
393 float right =
matrixEnergy( cluster, recHits, topology,
id, 1, 1, -2, 2 );
395 float centre =
matrixEnergy( cluster, recHits, topology,
id, 0, 0, -2, 2 );
398 return left > right ? left+centre : right+centre;
403 DetId id =
getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).
first;
406 float left =
matrixEnergy( cluster, recHits, topology,
id, -1, -1, -2, 2,flagsexcl, severitiesexcl, sevLv );
408 float right =
matrixEnergy( cluster, recHits, topology,
id, 1, 1, -2, 2,flagsexcl, severitiesexcl, sevLv );
410 float centre =
matrixEnergy( cluster, recHits, topology,
id, 0, 0, -2, 2,flagsexcl, severitiesexcl, sevLv );
413 return left > right ? left+centre : right+centre;
420 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, -2, 2 );
426 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, -2, 2 ,
427 flagsexcl, severitiesexcl, sevLv);
435 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, 0, 0 );
441 return matrixEnergy( cluster, recHits, topology,
id, -2, 2, 0, 0,flagsexcl, severitiesexcl, sevLv );
448 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, -1, 1 );
454 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, -1, 1,flagsexcl, severitiesexcl, sevLv );
461 return matrixEnergy( cluster, recHits, topology,
id, -1, 1, 0, 0 );
467 return matrixEnergy( cluster, recHits, topology,
id, -1, 1, 0, 0,flagsexcl, severitiesexcl, sevLv );
474 return matrixEnergy( cluster, recHits, topology,
id, -1, -1, 0, 0 );
480 return matrixEnergy( cluster, recHits, topology,
id, -1, -1, 0, 0,flagsexcl, severitiesexcl, sevLv );
487 return matrixEnergy( cluster, recHits, topology,
id, 1, 1, 0, 0 );
493 return matrixEnergy( cluster, recHits, topology,
id, 1, 1, 0, 0,flagsexcl, severitiesexcl, sevLv );
500 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, 1, 1 );
506 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, 1, 1,flagsexcl, severitiesexcl, sevLv );
514 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, -1, -1 );
520 return matrixEnergy( cluster, recHits, topology,
id, 0, 0, -1, -1,flagsexcl, severitiesexcl, sevLv );
526 float clusterEnergy = cluster.energy();
527 std::vector< std::pair<DetId, float> > v_id = cluster.hitsAndFractions();
529 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.";
530 return basketFraction;
532 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
535 std::sort( basketFraction.rbegin(), basketFraction.rend() );
536 return basketFraction;
543 float clusterEnergy = cluster.energy();
544 std::vector< std::pair<DetId, float> > v_id = cluster.hitsAndFractions();
546 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.";
547 return basketFraction;
549 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
552 std::sort( basketFraction.rbegin(), basketFraction.rend() );
553 return basketFraction;
559 float clusterEnergy = cluster.energy();
560 std::vector< std::pair<DetId, float> > v_id = cluster.hitsAndFractions();
562 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.";
563 return basketFraction;
565 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
568 std::sort( basketFraction.rbegin(), basketFraction.rend() );
569 return basketFraction;
576 float clusterEnergy = cluster.energy();
577 std::vector< std::pair<DetId, float> > v_id = cluster.hitsAndFractions();
579 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.";
580 return basketFraction;
582 for (
size_t i = 0;
i < v_id.size(); ++
i ) {
585 std::sort( basketFraction.rbegin(), basketFraction.rend() );
586 return basketFraction;
592 std::vector<EcalClusterTools::EcalClusterEnergyDeposition> energyDistribution;
595 CLHEP::Hep3Vector clVect(cluster.position().x(), cluster.position().y(), cluster.position().z());
596 CLHEP::Hep3Vector clDir(clVect);
597 clDir*=1.0/clDir.mag();
599 CLHEP::Hep3Vector theta_axis(clDir.y(),-clDir.x(),0.0);
600 theta_axis *= 1.0/theta_axis.mag();
601 CLHEP::Hep3Vector phi_axis = theta_axis.cross(clDir);
603 std::vector< std::pair<DetId, float> > clusterDetIds = cluster.hitsAndFractions();
607 std::vector< std::pair<DetId, float> >::iterator posCurrent;
609 for(posCurrent=clusterDetIds.begin(); posCurrent!=clusterDetIds.end(); ++posCurrent) {
613 if(( (*posCurrent).first !=
DetId(0)) && (recHits->
find( (*posCurrent).first ) != recHits->
end())) {
621 LogDebug(
"ClusterShapeAlgo") <<
"Crystal has insufficient energy: E = "
627 DetId id_ = (*posCurrent).first;
630 CLHEP::Hep3Vector gblPos (cellPos.
x(),cellPos.
y(),cellPos.
z());
632 CLHEP::Hep3Vector
diff = gblPos - clVect;
636 CLHEP::Hep3Vector DigiVect = diff - diff.dot(clDir)*clDir;
637 clEdep.
r = DigiVect.mag();
639 <<
"\tdiff = " << diff.mag()
640 <<
"\tr = " << clEdep.
r;
641 clEdep.
phi = DigiVect.angle(theta_axis);
642 if(DigiVect.dot(phi_axis)<0) clEdep.
phi = 2 *
M_PI - clEdep.
phi;
643 energyDistribution.push_back(clEdep);
646 return energyDistribution;
653 std::vector<EcalClusterTools::EcalClusterEnergyDeposition> energyDistribution =
getEnergyDepTopology( cluster, recHits, geometry, logW, w0 );
655 std::vector<float>
lat;
656 double r, redmoment=0;
657 double phiRedmoment = 0 ;
658 double etaRedmoment = 0 ;
660 int clusterSize=energyDistribution.size();
661 float etaLat_, phiLat_, lat_;
672 if (energyDistribution[1].deposited_energy >
673 energyDistribution[0].deposited_energy)
675 tmp=n2; n2=n1; n1=
tmp;
677 for (
int i=2;
i<clusterSize;
i++) {
679 if (energyDistribution[
i].deposited_energy >
680 energyDistribution[n1].deposited_energy)
683 n2 = n1; n1 =
i; n=
tmp;
685 if (energyDistribution[
i].deposited_energy >
686 energyDistribution[n2].deposited_energy)
692 r = energyDistribution[
n].r;
693 redmoment += r*r* energyDistribution[
n].deposited_energy;
694 double rphi = r *
cos (energyDistribution[n].
phi) ;
695 phiRedmoment += rphi * rphi * energyDistribution[
n].deposited_energy;
696 double reta = r *
sin (energyDistribution[n].phi) ;
697 etaRedmoment += reta * reta * energyDistribution[
n].deposited_energy;
699 double e1 = energyDistribution[n1].deposited_energy;
700 double e2 = energyDistribution[n2].deposited_energy;
702 lat_ = redmoment/(redmoment+2.19*2.19*(e1+e2));
703 phiLat_ = phiRedmoment/(phiRedmoment+2.19*2.19*(e1+e2));
704 etaLat_ = etaRedmoment/(etaRedmoment+2.19*2.19*(e1+e2));
706 lat.push_back(etaLat_);
707 lat.push_back(phiLat_);
719 for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
724 return meanPosition /
e5x5( cluster, recHits, topology );
734 for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
737 meanPosition = meanPosition +
recHitEnergy( *it, recHits,flagsexcl, severitiesexcl, sevLv ) *
position;
739 return meanPosition /
e5x5( cluster, recHits, topology,flagsexcl, severitiesexcl, sevLv );
757 std::vector<DetId> v_id =
matrixDetId( topology,seedId, -2, 2, -2, 2 );
758 for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
760 if(energy<0.)
continue;
767 return std::pair<float,float>(meanDEta,meanDPhi);
777 std::vector<DetId> v_id =
matrixDetId( topology,seedId, -2, 2, -2, 2 );
778 for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
780 if(energy<0.)
continue;
787 return std::pair<float,float>(meanDEta,meanDPhi);
800 std::pair<float,float> meanXY(0.,0.);
805 std::vector<DetId> v_id =
matrixDetId( topology,seedId, -2, 2, -2, 2 );
806 for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
808 if(energy<0.)
continue;
822 std::pair<float,float> meanXY(0.,0.);
827 std::vector<DetId> v_id =
matrixDetId( topology,seedId, -2, 2, -2, 2 );
828 for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
830 if(energy<0.)
continue;
843 float e_5x5 =
e5x5( cluster, recHits, topology );
844 float covEtaEta, covEtaPhi, covPhiPhi;
847 std::vector< std::pair<DetId, float> > v_id = cluster.hitsAndFractions();
851 double numeratorEtaEta = 0;
852 double numeratorEtaPhi = 0;
853 double numeratorPhiPhi = 0;
858 for (
int i = -2;
i <= 2; ++
i ) {
859 for (
int j = -2;
j <= 2; ++
j ) {
864 if ( energy <= 0 )
continue;
868 double dPhi = position.
phi() - meanPosition.phi();
872 double dEta = position.
eta() - meanPosition.eta();
877 numeratorEtaEta += w * dEta * dEta;
878 numeratorEtaPhi += w * dEta *
dPhi;
879 numeratorPhiPhi += w * dPhi *
dPhi;
883 if (denominator != 0.0) {
900 std::vector<float>
v;
901 v.push_back( covEtaEta );
902 v.push_back( covEtaPhi );
903 v.push_back( covPhiPhi );
911 float e_5x5 =
e5x5( cluster, recHits, topology,flagsexcl, severitiesexcl, sevLv );
912 float covEtaEta, covEtaPhi, covPhiPhi;
915 std::vector< std::pair<DetId, float> > v_id = cluster.hitsAndFractions();
919 double numeratorEtaEta = 0;
920 double numeratorEtaPhi = 0;
921 double numeratorPhiPhi = 0;
926 for (
int i = -2;
i <= 2; ++
i ) {
927 for (
int j = -2;
j <= 2; ++
j ) {
932 if ( energy <= 0 )
continue;
936 double dPhi = position.
phi() - meanPosition.phi();
940 double dEta = position.
eta() - meanPosition.eta();
945 numeratorEtaEta += w * dEta * dEta;
946 numeratorEtaPhi += w * dEta *
dPhi;
947 numeratorPhiPhi += w * dPhi *
dPhi;
951 if (denominator != 0.0) {
968 std::vector<float>
v;
969 v.push_back( covEtaEta );
970 v.push_back( covEtaPhi );
971 v.push_back( covPhiPhi );
986 float e_5x5 =
e5x5( cluster, recHits, topology );
987 float covEtaEta, covEtaPhi, covPhiPhi;
991 std::vector< std::pair<DetId, float> > v_id = cluster.hitsAndFractions();
996 double numeratorEtaEta = 0;
997 double numeratorEtaPhi = 0;
998 double numeratorPhiPhi = 0;
1003 const double barrelCrysSize = 0.01745;
1004 const double endcapCrysSize = 0.0447;
1009 const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1013 for (
int eastNr = -2; eastNr <= 2; ++eastNr ) {
1014 for (
int northNr = -2; northNr <= 2; ++northNr ) {
1018 if ( energy <= 0 )
continue;
1023 if(isBarrel) dPhi =
getNrCrysDiffInPhi(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.second;
1024 else dPhi =
getDPhiEndcap(*cursor,mean5x5XYPos.first,mean5x5XYPos.second);
1030 numeratorEtaEta += w * dEta * dEta;
1031 numeratorEtaPhi += w * dEta *
dPhi;
1032 numeratorPhiPhi += w * dPhi *
dPhi;
1038 if (denominator != 0.0) {
1039 covEtaEta = crysSize*crysSize* numeratorEtaEta /
denominator;
1040 covEtaPhi = crysSize*crysSize* numeratorEtaPhi /
denominator;
1041 covPhiPhi = crysSize*crysSize* numeratorPhiPhi /
denominator;
1056 std::vector<float>
v;
1057 v.push_back( covEtaEta );
1058 v.push_back( covEtaPhi );
1059 v.push_back( covPhiPhi );
1068 float e_5x5 =
e5x5( cluster, recHits, topology,flagsexcl, severitiesexcl, sevLv );
1069 float covEtaEta, covEtaPhi, covPhiPhi;
1073 std::vector< std::pair<DetId, float> > v_id = cluster.hitsAndFractions();
1075 std::pair<float,float> mean5x5XYPos =
mean5x5PositionInXY(cluster,recHits,topology,flagsexcl, severitiesexcl, sevLv);
1078 double numeratorEtaEta = 0;
1079 double numeratorEtaPhi = 0;
1080 double numeratorPhiPhi = 0;
1085 const double barrelCrysSize = 0.01745;
1086 const double endcapCrysSize = 0.0447;
1091 const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1095 for (
int eastNr = -2; eastNr <= 2; ++eastNr ) {
1096 for (
int northNr = -2; northNr <= 2; ++northNr ) {
1100 if ( energy <= 0 )
continue;
1105 if(isBarrel) dPhi =
getNrCrysDiffInPhi(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.second;
1106 else dPhi =
getDPhiEndcap(*cursor,mean5x5XYPos.first,mean5x5XYPos.second);
1112 numeratorEtaEta += w * dEta * dEta;
1113 numeratorEtaPhi += w * dEta *
dPhi;
1114 numeratorPhiPhi += w * dPhi *
dPhi;
1120 if (denominator != 0.0) {
1121 covEtaEta = crysSize*crysSize* numeratorEtaEta /
denominator;
1122 covEtaPhi = crysSize*crysSize* numeratorEtaPhi /
denominator;
1123 covPhiPhi = crysSize*crysSize* numeratorPhiPhi /
denominator;
1138 std::vector<float>
v;
1139 v.push_back( covEtaEta );
1140 v.push_back( covEtaPhi );
1141 v.push_back( covPhiPhi );
1163 if ((m>n) || ((n-m)%2 != 0) || (n<0) || (m<0))
return -1;
1167 if ((n>20) || (R0<=2.19))
return -1;
1175 double r,ph,
e,Re=0,Im=0;
1176 double TotalEnergy = cluster.energy();
1177 int index = (n/2)*(n/2)+(n/2)+m;
1178 std::vector<EcalClusterEnergyDeposition> energyDistribution =
getEnergyDepTopology( cluster, recHits, geometry, logW, w0 );
1179 int clusterSize = energyDistribution.size();
1180 if(clusterSize < 3)
return 0.0;
1182 for (
int i=0;
i<clusterSize;
i++)
1184 r = energyDistribution[
i].r / R0;
1186 std::vector<double> pol;
1187 pol.push_back(
f00(r) );
1188 pol.push_back(
f11(r) );
1189 pol.push_back(
f20(r) );
1190 pol.push_back(
f22(r) );
1191 pol.push_back(
f31(r) );
1192 pol.push_back(
f33(r) );
1193 pol.push_back(
f40(r) );
1194 pol.push_back(
f42(r) );
1195 pol.push_back(
f44(r) );
1196 pol.push_back(
f51(r) );
1197 pol.push_back(
f53(r) );
1198 pol.push_back(
f55(r) );
1199 ph = (energyDistribution[
i]).
phi;
1200 e = energyDistribution[
i].deposited_energy;
1201 Re = Re + e/TotalEnergy * pol[
index] *
cos( (
double) m * ph);
1202 Im = Im - e/TotalEnergy * pol[
index] *
sin( (
double) m * ph);
1205 return sqrt(Re*Re+Im*Im);
1212 double r, ph,
e, Re=0, Im=0, f_nm;
1213 double TotalEnergy = cluster.energy();
1214 std::vector<EcalClusterEnergyDeposition> energyDistribution =
getEnergyDepTopology( cluster, recHits, geometry, logW, w0 );
1215 int clusterSize=energyDistribution.size();
1216 if(clusterSize<3)
return 0.0;
1218 for (
int i = 0;
i < clusterSize; ++
i)
1220 r = energyDistribution[
i].r / R0;
1222 ph = energyDistribution[
i].phi;
1223 e = energyDistribution[
i].deposited_energy;
1225 for (
int s=0;
s<=(n-
m)/2;
s++) {
1232 Re = Re + e/TotalEnergy * f_nm *
cos( (
double) m*ph);
1233 Im = Im - e/TotalEnergy * f_nm *
sin( (
double) m*ph);
1236 return sqrt(Re*Re+Im*Im);
1253 return std::sqrt(iXNorm*iXNorm+iYNorm*iYNorm);
1280 int iXNorm = eeId.
ix()-50;
1281 if(iXNorm<=0) iXNorm--;
1292 int iYNorm = eeId.
iy()-50;
1293 if(iYNorm<=0) iYNorm--;
1302 float crysIEta =
getIEta(crysId);
1303 float orginIEta =
getIEta(orginId);
1306 float nrCrysDiff = crysIEta-orginIEta;
1311 if(crysIEta*orginIEta<0){
1312 if(crysIEta>0) nrCrysDiff--;
1322 float crysIPhi =
getIPhi(crysId);
1323 float orginIPhi =
getIPhi(orginId);
1326 float nrCrysDiff = crysIPhi-orginIPhi;
1329 if (nrCrysDiff > + 180) { nrCrysDiff = nrCrysDiff - 360; }
1330 if (nrCrysDiff < - 180) { nrCrysDiff = nrCrysDiff + 360; }
1341 float hitLocalR2 = (iXNorm-meanX)*(iXNorm-meanX)+(iYNorm-meanY)*(iYNorm-meanY);
1342 float hitR2 = iXNorm*iXNorm+iYNorm*iYNorm;
1343 float meanR2 = meanX*meanX+meanY*meanY;
1344 float hitR =
sqrt(hitR2);
1345 float meanR =
sqrt(meanR2);
1347 float tmp = (hitR2+meanR2-hitLocalR2)/(2*hitR*meanR);
1348 if (tmp<-1) tmp =-1;
1350 float phi = acos(tmp);
1360 float e_5x5 =
e5x5(bcluster, recHits, topology);
1361 float covEtaEta, covEtaPhi, covPhiPhi;
1368 double numeratorEtaEta = 0;
1369 double numeratorEtaPhi = 0;
1370 double numeratorPhiPhi = 0;
1373 const double barrelCrysSize = 0.01745;
1374 const double endcapCrysSize = 0.0447;
1379 const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1381 for (
size_t i = 0;
i < v_id.size(); ++
i) {
1385 if (energy <= 0)
continue;
1389 if(isBarrel) dPhi =
getNrCrysDiffInPhi(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.second;
1390 else dPhi =
getDPhiEndcap(*cursor,mean5x5XYPos.first,mean5x5XYPos.second);
1398 numeratorEtaEta += w * dEta * dEta;
1399 numeratorEtaPhi += w * dEta *
dPhi;
1400 numeratorPhiPhi += w * dPhi *
dPhi;
1404 if (denominator != 0.0) {
1405 covEtaEta = crysSize*crysSize* numeratorEtaEta /
denominator;
1406 covEtaPhi = crysSize*crysSize* numeratorEtaPhi /
denominator;
1407 covPhiPhi = crysSize*crysSize* numeratorPhiPhi /
denominator;
1422 std::vector<float>
v;
1423 v.push_back( covEtaEta );
1424 v.push_back( covEtaPhi );
1425 v.push_back( covPhiPhi );
1437 float e_5x5 =
e5x5(bcluster, recHits, topology);
1438 float covEtaEta, covEtaPhi, covPhiPhi;
1443 std::pair<float,float> mean5x5XYPos =
mean5x5PositionInXY(cluster,recHits,topology,flagsexcl, severitiesexcl, sevLv);
1445 double numeratorEtaEta = 0;
1446 double numeratorEtaPhi = 0;
1447 double numeratorPhiPhi = 0;
1450 const double barrelCrysSize = 0.01745;
1451 const double endcapCrysSize = 0.0447;
1456 const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1458 for (
size_t i = 0;
i < v_id.size(); ++
i) {
1462 if (energy <= 0)
continue;
1466 if(isBarrel) dPhi =
getNrCrysDiffInPhi(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.second;
1467 else dPhi =
getDPhiEndcap(*cursor,mean5x5XYPos.first,mean5x5XYPos.second);
1475 numeratorEtaEta += w * dEta * dEta;
1476 numeratorEtaPhi += w * dEta *
dPhi;
1477 numeratorPhiPhi += w * dPhi *
dPhi;
1481 if (denominator != 0.0) {
1482 covEtaEta = crysSize*crysSize* numeratorEtaEta /
denominator;
1483 covEtaPhi = crysSize*crysSize* numeratorEtaPhi /
denominator;
1484 covPhiPhi = crysSize*crysSize* numeratorPhiPhi /
denominator;
1499 std::vector<float>
v;
1500 v.push_back( covEtaEta );
1501 v.push_back( covEtaPhi );
1502 v.push_back( covPhiPhi );
1517 returnMoments.
sMaj = -1.;
1518 returnMoments.
sMin = -1.;
1519 returnMoments.
alpha = 0.;
1524 std::vector<const EcalRecHit*> RH_ptrs;
1526 std::vector< std::pair<DetId, float> > myHitsPair = basicCluster.hitsAndFractions();
1527 std::vector<DetId> usedCrystals;
1528 for(
unsigned int i=0;
i< myHitsPair.size();
i++){
1529 usedCrystals.push_back(myHitsPair[
i].
first);
1532 for(
unsigned int i=0;
i<usedCrystals.size();
i++){
1535 RH_ptrs.push_back( &(*myRH) );
1542 return returnMoments;
1551 returnMoments.
sMaj = -1.;
1552 returnMoments.
sMin = -1.;
1553 returnMoments.
alpha = 0.;
1560 return returnMoments;
1567 double mid_eta(0),mid_phi(0),mid_x(0),mid_y(0);
1571 double max_phi=-10.;
1572 double min_phi=100.;
1575 std::vector<double> etaDetId;
1576 std::vector<double> phiDetId;
1577 std::vector<double> xDetId;
1578 std::vector<double> yDetId;
1579 std::vector<double> wiDetId;
1581 unsigned int nCry=0;
1586 for(std::vector<const EcalRecHit*>::const_iterator rh_ptr = RH_ptrs.begin(); rh_ptr != RH_ptrs.end(); rh_ptr++){
1589 double temp_eta(0),temp_phi(0),temp_x(0),temp_y(0);
1590 isBarrel = (*rh_ptr)->detid().subdetId()==
EcalBarrel;
1593 temp_eta = (
getIEta((*rh_ptr)->detid()) > 0. ?
getIEta((*rh_ptr)->detid()) + 84.5 :
getIEta((*rh_ptr)->detid()) + 85.5);
1594 temp_phi=
getIPhi((*rh_ptr)->detid()) - 0.5;
1597 temp_eta =
getIEta((*rh_ptr)->detid());
1602 double temp_ene=(*rh_ptr)->energy();
1604 double temp_wi=((useLogWeights) ?
1609 if(temp_phi>max_phi) max_phi=temp_phi;
1610 if(temp_phi<min_phi) min_phi=temp_phi;
1611 etaDetId.push_back(temp_eta);
1612 phiDetId.push_back(temp_phi);
1613 xDetId.push_back(temp_x);
1614 yDetId.push_back(temp_y);
1615 wiDetId.push_back(temp_wi);
1616 denominator+=temp_wi;
1622 if(max_phi==359.5 && min_phi==0.5){
1623 for(
unsigned int i=0;
i<nCry;
i++){
1624 if(phiDetId[
i] - 179. > 0.) phiDetId[
i]-=360.;
1625 mid_phi+=phiDetId[
i]*wiDetId[
i];
1626 mid_eta+=etaDetId[
i]*wiDetId[
i];
1629 for(
unsigned int i=0;
i<nCry;
i++){
1630 mid_phi+=phiDetId[
i]*wiDetId[
i];
1631 mid_eta+=etaDetId[
i]*wiDetId[
i];
1635 for(
unsigned int i=0;
i<nCry;
i++){
1636 mid_eta+=etaDetId[
i]*wiDetId[
i];
1637 mid_x+=xDetId[
i]*wiDetId[
i];
1638 mid_y+=yDetId[
i]*wiDetId[
i];
1654 double deta(0),dphi(0);
1656 for(
unsigned int i=0;
i<nCry;
i++) {
1658 deta = etaDetId[
i]-mid_eta;
1659 dphi = phiDetId[
i]-mid_phi;
1661 deta = etaDetId[
i]-mid_eta;
1662 float hitLocalR2 = (xDetId[
i]-mid_x)*(xDetId[
i]-mid_x)+(yDetId[
i]-mid_y)*(yDetId[
i]-mid_y);
1663 float hitR2 = xDetId[
i]*xDetId[
i]+yDetId[
i]*yDetId[
i];
1664 float meanR2 = mid_x*mid_x+mid_y*mid_y;
1665 float hitR =
sqrt(hitR2);
1666 float meanR =
sqrt(meanR2);
1667 float phi = acos((hitR2+meanR2-hitLocalR2)/(2*hitR*meanR));
1671 See += (wiDetId[
i]* deta * deta) / denominator;
1672 Spp += phiCorrectionFactor*(wiDetId[
i]* dphi * dphi) / denominator;
1673 Sep +=
sqrt(phiCorrectionFactor)*(wiDetId[
i]*deta*dphi) / denominator;
1679 returnMoments.
sMaj = ((See + Spp) +
sqrt((See - Spp)*(See - Spp) + 4.*Sep*Sep)) / 2.;
1680 returnMoments.
sMin = ((See + Spp) -
sqrt((See - Spp)*(See - Spp) + 4.*Sep*Sep)) / 2.;
1682 returnMoments.
alpha = atan( (See - Spp +
sqrt( (Spp - See)*(Spp - See) + 4.*Sep*Sep )) / (2.*Sep));
1684 return returnMoments;
1699 std::vector<const EcalRecHit*>RH_ptrs;
1700 std::vector< std::pair<DetId, float> > myHitsPair = superCluster.
hitsAndFractions();
1701 std::vector<DetId> usedCrystals;
1702 for(
unsigned int i=0;
i< myHitsPair.size();
i++){
1703 usedCrystals.push_back(myHitsPair[
i].
first);
1705 for(
unsigned int i=0;
i<usedCrystals.size();
i++){
1708 if( myRH != recHits.
end() && myRH->energy() > energyThreshold){
1709 RH_ptrs.push_back( &(*myRH) );
1723 std::vector<const EcalRecHit*>RH_ptrs;
1724 std::vector< std::pair<DetId, float> > myHitsPair = superCluster.
hitsAndFractions();
1725 std::vector<DetId> usedCrystals;
1726 for(
unsigned int i=0;
i< myHitsPair.size();
i++){
1727 usedCrystals.push_back(myHitsPair[
i].
first);
1730 for(
unsigned int i=0;
i<usedCrystals.size();
i++){
1733 if(myRH != recHits.
end() && myRH->energy() > energyRHThresh)
1734 RH_ptrs.push_back( &(*myRH) );
1741 EBDetId EBdetIdi( rh->detid() );
1743 bool inEtaWindow = (
abs(
deltaIEta(seedPosition[0],EBdetIdi.ieta()) ) <= ieta_delta );
1744 bool inPhiWindow = (
abs(
deltaIPhi(seedPosition[1],EBdetIdi.iphi()) ) <= iphi_delta );
1745 bool passEThresh = ( rh->energy() > energyRHThresh );
1746 bool alreadyCounted =
false;
1749 bool is_SCrh_inside_recHits =
false;
1750 for(
unsigned int i=0;
i<usedCrystals.size();
i++){
1752 if(SCrh != recHits.
end()){
1753 is_SCrh_inside_recHits =
true;
1754 if( rh->detid() == SCrh->detid() ) alreadyCounted =
true;
1758 if( is_SCrh_inside_recHits && !alreadyCounted && passEThresh && inEtaWindow && inPhiWindow){
1759 RH_ptrs.push_back( &(*rh) );
1772 std::vector<float> shapes;
1775 if(RH_ptrs.size()<2){
1776 shapes.push_back( -3 );
1777 shapes.push_back( -3 );
1783 int tempInt = seedPosition[0];
1784 if(tempInt <0) tempInt++;
1788 float centerIEta = 0.;
1789 float centerIPhi = 0.;
1792 for(std::vector<const EcalRecHit*>::const_iterator rh_ptr = RH_ptrs.begin(); rh_ptr != RH_ptrs.end(); rh_ptr++){
1794 EBDetId EBdetIdi( (*rh_ptr)->detid() );
1795 if(fabs(energyTotal) < 0.0001){
1797 shapes.push_back( -2 );
1798 shapes.push_back( -2 );
1802 if(fabs(weightedPositionMethod)<0.0001){
1803 weight = (*rh_ptr)->energy()/energyTotal;
1805 weight =
std::max(0.0, 4.2 +
log((*rh_ptr)->energy()/energyTotal));
1808 centerIEta += weight*
deltaIEta(seedPosition[0],EBdetIdi.ieta());
1809 centerIPhi += weight*
deltaIPhi(seedPosition[1],EBdetIdi.iphi());
1811 if(fabs(denominator) < 0.0001){
1813 shapes.push_back( -2 );
1814 shapes.push_back( -2 );
1822 TMatrixDSym inertia(2);
1823 double inertia00 = 0.;
1824 double inertia01 = 0.;
1825 double inertia11 = 0.;
1827 for(std::vector<const EcalRecHit*>::const_iterator rh_ptr = RH_ptrs.begin(); rh_ptr != RH_ptrs.end(); rh_ptr++){
1829 EBDetId EBdetIdi( (*rh_ptr)->detid() );
1831 if(fabs(energyTotal) < 0.0001){
1833 shapes.push_back( -2 );
1834 shapes.push_back( -2 );
1838 if(fabs(weightedPositionMethod) < 0.0001){
1839 weight = (*rh_ptr)->energy()/energyTotal;
1841 weight =
std::max(0.0, 4.2 +
log((*rh_ptr)->energy()/energyTotal));
1844 float ieta_rh_to_center =
deltaIEta(seedPosition[0],EBdetIdi.ieta()) - centerIEta;
1845 float iphi_rh_to_center =
deltaIPhi(seedPosition[1],EBdetIdi.iphi()) - centerIPhi;
1847 inertia00 += weight*iphi_rh_to_center*iphi_rh_to_center;
1848 inertia01 -= weight*iphi_rh_to_center*ieta_rh_to_center;
1849 inertia11 += weight*ieta_rh_to_center*ieta_rh_to_center;
1853 inertia[0][0] = inertia00;
1854 inertia[0][1] = inertia01;
1855 inertia[1][0] = inertia01;
1856 inertia[1][1] = inertia11;
1860 TMatrixD eVectors(2,2);
1861 TVectorD eValues(2);
1863 eVectors=inertia.EigenVectors(eValues);
1870 TVectorD smallerAxis(2);
1871 smallerAxis[0]=eVectors[0][1];
1872 smallerAxis[1]=eVectors[1][1];
1875 Double_t
temp = fabs(smallerAxis[0]);
1876 if(fabs(eValues[0]) < 0.0001){
1878 shapes.push_back( -2 );
1879 shapes.push_back( -2 );
1883 float Roundness = eValues[1]/eValues[0];
1884 float Angle=acos(temp);
1886 if( -0.00001 < Roundness && Roundness < 0) Roundness = 0.;
1887 if( -0.00001 < Angle && Angle < 0 ) Angle = 0.;
1889 shapes.push_back( Roundness );
1890 shapes.push_back( Angle );
1899 int rel_iphi = rh_iphi - seed_iphi;
1901 if(rel_iphi > 180) rel_iphi = rel_iphi - 360;
1902 if(rel_iphi < -180) rel_iphi = rel_iphi + 360;
1911 if(seed_ieta < 0) seed_ieta++;
1912 if(rh_ieta < 0) rh_ieta++;
1913 int rel_ieta = rh_ieta - seed_ieta;
1919 std::vector<int> seedPosition;
1924 for(std::vector<const EcalRecHit*>::const_iterator rh_ptr = RH_ptrs.begin(); rh_ptr != RH_ptrs.end(); rh_ptr++){
1927 EBDetId EBdetIdi( (*rh_ptr)->detid() );
1929 if(eSeedRH < (*rh_ptr)->energy()){
1930 eSeedRH = (*rh_ptr)->energy();
1931 iEtaSeedRH = EBdetIdi.
ieta();
1932 iPhiSeedRH = EBdetIdi.iphi();
1937 seedPosition.push_back(iEtaSeedRH);
1938 seedPosition.push_back(iPhiSeedRH);
1939 return seedPosition;
1947 for(std::vector<const EcalRecHit*>::const_iterator rh_ptr = RH_ptrs.begin(); rh_ptr != RH_ptrs.end(); rh_ptr++){
1948 sumE += (*rh_ptr)->energy();
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
void home() const
move the navigator back to the starting point
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
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
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)
virtual T offsetBy(int deltaX, int deltaY) const
Free movement of arbitray steps.
Cos< T >::type cos(const T &t)
static const int kModulesPerSM
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)
static int position[264][3]
const CaloClusterPtr & seed() const
seed BasicCluster
double energySum(const DataFrame &df, int fs, int ls)
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