EcalSeverityLevelAlgo Class Reference

#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 &)

Detailed Description

Definition at line 8 of file EcalSeverityLevelAlgo.h.

Member Enumeration Documentation

enum EcalSeverityLevelAlgo::EcalSeverityLevel

Enumerator
kGood
kProblematic
kRecovered
kTime
kWeird
kBad

Definition at line 19 of file EcalSeverityLevelAlgo.h.

{ kGood=0, kProblematic, kRecovered, kTime, kWeird, kBad };

enum EcalSeverityLevelAlgo::SpikeId

Enumerator
kE1OverE9
kSwissCross
kSwissCrossBordersIncluded

Definition at line 21 of file EcalSeverityLevelAlgo.h.

{ kE1OverE9=0, kSwissCross, kSwissCrossBordersIncluded };

Member Function Documentation

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  )

staticprivate

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  )

staticprivate

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  )

staticprivate

Definition at line 201 of file EcalSeverityLevelAlgo.cc.

References EcalBarrel, EcalEndcap, EBDetId::offsetBy(), EEDetId::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  )

staticprivate

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 ( 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 EBOccupancyTask::analyze(), EBTimingTask::analyze(), EETimingTask::analyze(), EEOccupancyTask::analyze(), QcdPhotonsDQM::analyze(), EcalClusterSeverityLevelAlgo::closestProblematic(), egammaisolation::EgammaRecHitExtractor::collect(), CaloTowersCreationAlgo::convert(), CaloTowersCreationAlgo::ecalChanStatusForCaloTower(), spr::eECALmatrix(), EgammaRecHitIsolation::getSum_(), EcalClusterSeverityLevelAlgo::goodFraction(), HybridClusterAlgo::makeClusters(), EleIsoDetIdCollectionProducer::produce(), GamIsoDetIdCollectionProducer::produce(), HiSpikeCleaner::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;
}

int EcalSeverityLevelAlgo::severityLevel ( uint32_t  rhFlag, uint16_t  dbStatus  )

staticprivate

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  )

staticprivate

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 );
}

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 EgammaHLTR9Producer::produce(), HiSpikeCleaner::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;
}