|
|
|
#include <EcalSeverityLevelAlgo.h>
Public Types | |
| enum | EcalSeverityLevel { kGood = 0, kProblematic, kRecovered, kTime, kWeird, kBad } |
| enum | SpikeId { kE1OverE9 = 0, kSwissCross, kSwissCrossBordersIncluded } |
Static Public Member Functions | |
| static float | E1OverE9 (const DetId id, const EcalRecHitCollection &, float recHitEtThreshold=0.) |
| static float | E2overE9 (const DetId id, const EcalRecHitCollection &, float recHitEtThreshold=10.0, float recHitEtThreshold2=1.0, bool avoidIeta85=false, bool KillSecondHit=true) |
| static int | severityLevel (const DetId, const EcalRecHitCollection &, const EcalChannelStatus &, float recHitEtThreshold=5., SpikeId spId=kSwissCross, float spIdThreshold=0.95, float recHitEnergyThresholdForTiming=2., float recHitEnergyThresholdForEE=15, float spIdThresholdIEta85=0.999) |
| static float | spikeFromNeighbours (const DetId id, const EcalRecHitCollection &, float recHitEtThreshold, SpikeId spId) |
| static bool | spikeFromTiming (const EcalRecHit &, float recHitEnergyThreshold) |
| static float | swissCross (const DetId id, const EcalRecHitCollection &, float recHitEtThreshold=0., bool avoidIeta85=true) |
Static Private Member Functions | |
| static float | recHitApproxEt (const DetId id, const EcalRecHitCollection &recHits) |
| static float | recHitE (const DetId id, const EcalRecHitCollection &recHits) |
| static float | recHitE (const DetId id, const EcalRecHitCollection &recHits, int dEta, int dPhi) |
| static uint16_t | retrieveDBStatus (const DetId, const EcalChannelStatus &chStatus) |
| static int | severityLevel (uint32_t rhFlag, uint16_t dbStatus) |
| static int | severityLevel (const EcalRecHit &, const EcalChannelStatus &) |
Definition at line 8 of file EcalSeverityLevelAlgo.h.
Definition at line 19 of file EcalSeverityLevelAlgo.h.
{ kGood=0, kProblematic, kRecovered, kTime, kWeird, kBad };
Definition at line 21 of file EcalSeverityLevelAlgo.h.
| float EcalSeverityLevelAlgo::E1OverE9 | ( | const DetId | id, |
| const EcalRecHitCollection & | recHits, | ||
| float | recHitEtThreshold = 0. |
||
| ) | [static] |
ratio between the crystal energy and the energy in the 3x3 matrix of crystal
Definition at line 133 of file EcalSeverityLevelAlgo.cc.
References EcalBarrel, EcalEndcap, recHitApproxEt(), and recHitE().
Referenced by spikeFromNeighbours().
{
// compute E1 over E9
if ( id.subdetId() == EcalBarrel ) {
// select recHits with Et above recHitEtThreshold
if ( recHitApproxEt( id, recHits ) < recHitEtThreshold ) return 0;
EBDetId ebId( id );
float s9 = 0;
for ( int deta = -1; deta <= +1; ++deta ) {
for ( int dphi = -1; dphi <= +1; ++dphi ) {
s9 += recHitE( id, recHits, deta, dphi );
}
}
return recHitE(id, recHits) / s9;
} else if( id.subdetId() == EcalEndcap ) {
// select recHits with Et above recHitEtThreshold
if ( recHitApproxEt( id, recHits ) < recHitEtThreshold ) return 0;
EEDetId eeId( id );
float s9 = 0;
for ( int dx = -1; dx <= +1; ++dx ) {
for ( int dy = -1; dy <= +1; ++dy ) {
s9 += recHitE( id, recHits, dx, dy );
}
}
return recHitE(id, recHits) / s9;
}
return 0;
}
| float EcalSeverityLevelAlgo::E2overE9 | ( | const DetId | id, |
| const EcalRecHitCollection & | recHits, | ||
| float | recHitEtThreshold = 10.0, |
||
| float | recHitEtThreshold2 = 1.0, |
||
| bool | avoidIeta85 = false, |
||
| bool | KillSecondHit = true |
||
| ) | [static] |
| | | | +-+-+-+ | |1|2| +-+-+-+ | | | |
1 - input hit, 2 - highest energy hit in a 3x3 around 1
rechit 1 must have E_t > recHitEtThreshold rechit 2 must have E_t > recHitEtThreshold2
function returns value of E2/E9 centered around 1 (E2=energy of hits 1+2) if energy of 1>2
if energy of 2>1 and KillSecondHit is set to true, function returns value of E2/E9 centered around 2 *provided* that 1 is the highest energy hit in a 3x3 centered around 2, otherwise, function returns 0
Definition at line 245 of file EcalSeverityLevelAlgo.cc.
References abs, EcalBarrel, EcalEndcap, EBDetId::ieta(), min, EBDetId::offsetBy(), recHitApproxEt(), and recHitE().
{
// compute e2overe9
//
// | | | |
// +-+-+-+
// | |1|2|
// +-+-+-+
// | | | |
//
// 1 - input hit, 2 - highest energy hit in a 3x3 around 1
//
// rechit 1 must have E_t > recHitEtThreshold
// rechit 2 must have E_t > recHitEtThreshold2
//
// function returns value of E2/E9 centered around 1 (E2=energy of hits 1+2) if energy of 1>2
//
// if energy of 2>1 and KillSecondHit is set to true, function returns value of E2/E9 centered around 2
// *provided* that 1 is the highest energy hit in a 3x3 centered around 2, otherwise, function returns 0
if ( id.subdetId() == EcalBarrel ) {
EBDetId ebId( id );
// avoid recHits at |eta|=85 where one side of the neighbours is missing
if ( abs(ebId.ieta())==85 && avoidIeta85) return 0;
// select recHits with Et above recHitEtThreshold
float e1 = recHitE( id, recHits );
float ete1=recHitApproxEt( id, recHits );
// check that rechit E_t is above threshold
if (ete1 < std::min(recHitEtThreshold,recHitEtThreshold2) ) return 0;
if (ete1 < recHitEtThreshold && !KillSecondHit ) return 0;
float e2=-1;
float ete2=0;
float s9 = 0;
// coordinates of 2nd hit relative to central hit
int e2eta=0;
int e2phi=0;
// LOOP OVER 3x3 ARRAY CENTERED AROUND HIT 1
for ( int deta = -1; deta <= +1; ++deta ) {
for ( int dphi = -1; dphi <= +1; ++dphi ) {
// compute 3x3 energy
float etmp=recHitE( id, recHits, deta, dphi );
s9 += etmp;
EBDetId idtmp=EBDetId::offsetBy(id,deta,dphi);
float eapproxet=recHitApproxEt( idtmp, recHits );
// remember 2nd highest energy deposit (above threshold) in 3x3 array
if (etmp>e2 && eapproxet>recHitEtThreshold2 && !(deta==0 && dphi==0)) {
e2=etmp;
ete2=eapproxet;
e2eta=deta;
e2phi=dphi;
}
}
}
if ( e1 == 0 ) return 0;
// return 0 if 2nd hit is below threshold
if ( e2 == -1 ) return 0;
// compute e2/e9 centered around 1st hit
float e2nd=e1+e2;
float e2e9=0;
if (s9!=0) e2e9=e2nd/s9;
// if central hit has higher energy than 2nd hit
// return e2/e9 if 1st hit is above E_t threshold
if (e1 > e2 && ete1>recHitEtThreshold) return e2e9;
// if second hit has higher energy than 1st hit
if ( e2 > e1 ) {
// return 0 if user does not want to flag 2nd hit, or
// hits are below E_t thresholds - note here we
// now assume the 2nd hit to be the leading hit.
if (!KillSecondHit || ete2<recHitEtThreshold || ete1<recHitEtThreshold2) {
return 0;
}
else {
// LOOP OVER 3x3 ARRAY CENTERED AROUND HIT 2
float s92nd=0;
float e2nd_prime=0;
int e2prime_eta=0;
int e2prime_phi=0;
EBDetId secondid=EBDetId::offsetBy(id,e2eta,e2phi);
for ( int deta = -1; deta <= +1; ++deta ) {
for ( int dphi = -1; dphi <= +1; ++dphi ) {
// compute 3x3 energy
float etmp=recHitE( secondid, recHits, deta, dphi );
s92nd += etmp;
if (etmp>e2nd_prime && !(deta==0 && dphi==0)) {
e2nd_prime=etmp;
e2prime_eta=deta;
e2prime_phi=dphi;
}
}
}
// if highest energy hit around E2 is not the same as the input hit, return 0;
if (!(e2prime_eta==-e2eta && e2prime_phi==-e2phi))
{
return 0;
}
// compute E2/E9 around second hit
float e2e9_2=0;
if (s92nd!=0) e2e9_2=e2nd/s92nd;
// return the value of E2/E9 calculated around 2nd hit
return e2e9_2;
}
}
} else if ( id.subdetId() == EcalEndcap ) {
// only used for EB at the moment
return 0;
}
return 0;
}
| float EcalSeverityLevelAlgo::recHitApproxEt | ( | const DetId | id, |
| const EcalRecHitCollection & | recHits | ||
| ) | [static, private] |
Definition at line 226 of file EcalSeverityLevelAlgo.cc.
References EBDetId::approxEta(), EcalBarrel, and recHitE().
Referenced by E1OverE9(), E2overE9(), and swissCross().
{
// for the time being works only for the barrel
if ( id.subdetId() == EcalBarrel ) {
return recHitE( id, recHits ) / cosh( EBDetId::approxEta( id ) );
}
return 0;
}
| float EcalSeverityLevelAlgo::recHitE | ( | const DetId | id, |
| const EcalRecHitCollection & | recHits | ||
| ) | [static, private] |
return energy of a recHit (if in the collection)
Definition at line 214 of file EcalSeverityLevelAlgo.cc.
References edm::SortedCollection< T, SORT >::end(), and edm::SortedCollection< T, SORT >::find().
Referenced by E1OverE9(), E2overE9(), recHitApproxEt(), recHitE(), severityLevel(), and swissCross().
{
if ( id == DetId(0) ) {
return 0;
} else {
EcalRecHitCollection::const_iterator it = recHits.find( id );
if ( it != recHits.end() ) return (*it).energy();
}
return 0;
}
| float EcalSeverityLevelAlgo::recHitE | ( | const DetId | id, |
| const EcalRecHitCollection & | recHits, | ||
| int | dEta, | ||
| int | dPhi | ||
| ) | [static, private] |
Definition at line 201 of file EcalSeverityLevelAlgo.cc.
References EcalBarrel, EcalEndcap, EEDetId::offsetBy(), EBDetId::offsetBy(), and recHitE().
{
// in the barrel: di = dEta dj = dPhi
// in the endcap: di = dX dj = dY
DetId nid;
if( id.subdetId() == EcalBarrel) nid = EBDetId::offsetBy( id, di, dj );
else if( id.subdetId() == EcalEndcap) nid = EEDetId::offsetBy( id, di, dj );
return ( nid == DetId(0) ? 0 : recHitE( nid, recHits ) );
}
| uint16_t EcalSeverityLevelAlgo::retrieveDBStatus | ( | const DetId | id, |
| const EcalChannelStatus & | chStatus | ||
| ) | [static, private] |
Definition at line 94 of file EcalSeverityLevelAlgo.cc.
References EcalCondObjectContainer< T >::end(), and EcalCondObjectContainer< T >::find().
Referenced by severityLevel().
{
EcalChannelStatus::const_iterator chIt = chStatus.find( id );
uint16_t dbStatus = 0;
if ( chIt != chStatus.end() ) {
dbStatus = chIt->getStatusCode();
} else {
edm::LogError("EcalSeverityLevelError") << "No channel status found for xtal "
<< id.rawId()
<< "! something wrong with EcalChannelStatus in your DB? ";
}
return dbStatus;
}
| int EcalSeverityLevelAlgo::severityLevel | ( | uint32_t | rhFlag, |
| uint16_t | dbStatus | ||
| ) | [static, private] |
Definition at line 64 of file EcalSeverityLevelAlgo.cc.
References kBad, EcalRecHit::kDead, EcalRecHit::kFakeNeighbours, EcalRecHit::kFaultyHardware, kGood, EcalRecHit::kKilled, EcalRecHit::kLeadingEdgeRecovered, EcalRecHit::kNeighboursRecovered, EcalRecHit::kNoisy, EcalRecHit::kOutOfTime, EcalRecHit::kPoorCalib, EcalRecHit::kPoorReco, kProblematic, kRecovered, EcalRecHit::kSaturated, and EcalRecHit::kTowerRecovered.
{
// DB info currently not used at this level
if ( rhFlag == EcalRecHit::kPoorReco
|| rhFlag == EcalRecHit::kOutOfTime
|| rhFlag == EcalRecHit::kNoisy
|| rhFlag == EcalRecHit::kPoorCalib
|| rhFlag == EcalRecHit::kFaultyHardware
) {
// problematic
return kProblematic;
} else if ( rhFlag == EcalRecHit::kLeadingEdgeRecovered
|| rhFlag == EcalRecHit::kNeighboursRecovered
|| rhFlag == EcalRecHit::kTowerRecovered
) {
// recovered
return kRecovered;
} else if ( rhFlag == EcalRecHit::kDead
|| rhFlag == EcalRecHit::kSaturated
//|| rhFlag == EcalRecHit::kFake // will be uncommented when validated
|| rhFlag == EcalRecHit::kFakeNeighbours
|| rhFlag == EcalRecHit::kKilled ) {
// recovery failed (or not tried) or signal is fake or channel
// is dead
return kBad;
}
// good
return kGood;
}
| int EcalSeverityLevelAlgo::severityLevel | ( | const EcalRecHit & | recHit, |
| const EcalChannelStatus & | chStatus | ||
| ) | [static, private] |
Definition at line 54 of file EcalSeverityLevelAlgo.cc.
References EcalRecHit::id(), EcalRecHit::recoFlag(), retrieveDBStatus(), and severityLevel().
{
// the channel is there, check its flags
// and combine with DB (not needed at the moment)
uint32_t rhFlag = recHit.recoFlag();
uint16_t dbStatus = retrieveDBStatus( recHit.id(), chStatus );
return severityLevel( rhFlag, dbStatus );
}
| int EcalSeverityLevelAlgo::severityLevel | ( | const DetId | id, |
| const EcalRecHitCollection & | recHits, | ||
| const EcalChannelStatus & | chStatus, | ||
| float | recHitEtThreshold = 5., |
||
| SpikeId | spId = kSwissCross, |
||
| float | spIdThreshold = 0.95, |
||
| float | recHitEnergyThresholdForTiming = 2., |
||
| float | recHitEnergyThresholdForEE = 15, |
||
| float | spIdThresholdIEta85 = 0.999 |
||
| ) | [static] |
compute the severity level
Definition at line 5 of file EcalSeverityLevelAlgo.cc.
References abs, EcalBarrel, EcalEndcap, edm::SortedCollection< T, SORT >::end(), edm::SortedCollection< T, SORT >::find(), kBad, kGood, kProblematic, kSwissCross, kSwissCrossBordersIncluded, kTime, kWeird, funct::log(), recHitE(), retrieveDBStatus(), spikeFromNeighbours(), spikeFromTiming(), and swissCross().
Referenced by EEOccupancyTask::analyze(), EETimingTask::analyze(), EBOccupancyTask::analyze(), QcdPhotonsDQM::analyze(), EBTimingTask::analyze(), EcalClusterSeverityLevelAlgo::closestProblematic(), egammaisolation::EgammaRecHitExtractor::collect(), CaloTowersCreationAlgo::convert(), CaloTowersCreationAlgo::ecalChanStatusForCaloTower(), spr::eECALmatrix(), EgammaRecHitIsolation::getSum_(), EcalClusterSeverityLevelAlgo::goodFraction(), HybridClusterAlgo::makeClusters(), GamIsoDetIdCollectionProducer::produce(), HiSpikeCleaner::produce(), EleIsoDetIdCollectionProducer::produce(), severityLevel(), and ObjectValidator::validHit().
{
// get DB flag
uint16_t dbStatus = retrieveDBStatus( id, chStatus );
// get recHit flags
EcalRecHitCollection::const_iterator it = recHits.find( id );
if ( it == recHits.end() ) {
// the channel is not in the recHit collection:
// dead or zero-suppressed?
if ( dbStatus >= 10 ) { // originally dead
return kBad;
} else if ( dbStatus > 0 && dbStatus < 10 ) {
// zero-suppressed and originally problematic
return kProblematic;
} else if ( dbStatus == 0 ) {
// zero-suppressed and originally good
return kGood;
}
} else {
// the channel is in the recHit collection
// .. is it a spike?
// check the topology
if ( id.subdetId() == EcalBarrel ) {
if ( abs(((EBDetId)id).ieta()) == 85 && spId == kSwissCrossBordersIncluded && spikeFromNeighbours(id, recHits, recHitEtThreshold, spId) > spIdThresholdIEta85 ) return kWeird;
if ( spikeFromNeighbours(id, recHits, recHitEtThreshold, spId) > spIdThreshold ) return kWeird;
}
// check the timing (currently only a trivial check)
if ( id.subdetId() == EcalBarrel && spikeFromTiming(*it, recHitEnergyThresholdForTiming) ) return kTime;
// filtering on VPT discharges in the endcap
// introduced >= kSwisscross to take borders too, SA 20100913
float re = 0; // protect the log function: make the computation possible only for re > 0
if ( id.subdetId() == EcalEndcap && spId >= kSwissCross && (re = recHitE(id, recHits)) > 0 && ( 1-swissCross(id, recHits, recHitEnergyThresholdForEE, spId) < 0.02*log(re/4.) ) ) return kWeird;
// .. not a spike, return the normal severity level
return severityLevel( *it, chStatus );
}
return kGood;
}
| float EcalSeverityLevelAlgo::spikeFromNeighbours | ( | const DetId | id, |
| const EcalRecHitCollection & | recHits, | ||
| float | recHitEtThreshold, | ||
| SpikeId | spId | ||
| ) | [static] |
return the estimator of the signal being a spike based on the topological information from the neighbours
Definition at line 108 of file EcalSeverityLevelAlgo.cc.
References E1OverE9(), kE1OverE9, kSwissCross, kSwissCrossBordersIncluded, and swissCross().
Referenced by severityLevel().
{
switch( spId ) {
case kE1OverE9:
return E1OverE9( id, recHits, recHitThreshold );
break;
case kSwissCross:
return swissCross( id, recHits, recHitThreshold , true);
break;
case kSwissCrossBordersIncluded:
return swissCross( id, recHits, recHitThreshold , false);
break;
default:
edm::LogInfo("EcalSeverityLevelAlgo") << "Algorithm number " << spId
<< " not known. Please check the enum in EcalSeverityLevelAlgo.h";
break;
}
return 0;
}
| bool EcalSeverityLevelAlgo::spikeFromTiming | ( | const EcalRecHit & | recHit, |
| float | recHitEnergyThreshold | ||
| ) | [static] |
return whether or not the rechit is a spike based on the kOutOfTime rechit flag
Definition at line 236 of file EcalSeverityLevelAlgo.cc.
References CaloRecHit::energy(), EcalRecHit::kOutOfTime, and EcalRecHit::recoFlag().
Referenced by severityLevel().
{
if ( recHit.energy() < recHitEnergyThreshold ) return false;
if ( recHit.recoFlag() == EcalRecHit::kOutOfTime ) return true;
return false;
}
| float EcalSeverityLevelAlgo::swissCross | ( | const DetId | id, |
| const EcalRecHitCollection & | recHits, | ||
| float | recHitEtThreshold = 0., |
||
| bool | avoidIeta85 = true |
||
| ) | [static] |
1 - the ratio between the energy in the swiss cross around a crystal and the crystal energy (also called S4/S1, Rook)
Definition at line 163 of file EcalSeverityLevelAlgo.cc.
References abs, EcalBarrel, EcalEndcap, EBDetId::ieta(), recHitApproxEt(), and recHitE().
Referenced by HiSpikeCleaner::produce(), EgammaHLTR9Producer::produce(), severityLevel(), and spikeFromNeighbours().
{
// compute swissCross
if ( id.subdetId() == EcalBarrel ) {
EBDetId ebId( id );
// avoid recHits at |eta|=85 where one side of the neighbours is missing
// (may improve considering also eta module borders, but no
// evidence for the time being that there the performance is
// different)
if ( abs(ebId.ieta())==85 && avoidIeta85) return 0;
// select recHits with Et above recHitThreshold
if ( recHitApproxEt( id, recHits ) < recHitThreshold ) return 0;
float s4 = 0;
float e1 = recHitE( id, recHits );
// protect against nan (if 0 threshold is given above)
if ( e1 == 0 ) return 0;
s4 += recHitE( id, recHits, 1, 0 );
s4 += recHitE( id, recHits, -1, 0 );
s4 += recHitE( id, recHits, 0, 1 );
s4 += recHitE( id, recHits, 0, -1 );
return 1 - s4 / e1;
} else if ( id.subdetId() == EcalEndcap ) {
EEDetId eeId( id );
// select recHits with E above recHitThreshold
float e1 = recHitE( id, recHits );
if ( e1 < recHitThreshold ) return 0;
float s4 = 0;
// protect against nan (if 0 threshold is given above)
if ( e1 == 0 ) return 0;
s4 += recHitE( id, recHits, 1, 0 );
s4 += recHitE( id, recHits, -1, 0 );
s4 += recHitE( id, recHits, 0, 1 );
s4 += recHitE( id, recHits, 0, -1 );
return 1 - s4 / e1;
}
return 0;
}
1.7.3