CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoEcal/EgammaCoreTools/src/EcalTools.cc

Go to the documentation of this file.
00001 // $Id: EcalTools.cc,v 1.4 2011/05/19 14:39:29 argiro Exp $
00002 
00003 #include "RecoEcal/EgammaCoreTools/interface/EcalTools.h"
00004 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00005 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00006 #include "RecoEcal/EgammaCoreTools/interface/EcalNextToDeadChannel.h"
00007 #include "RecoEcal/EgammaCoreTools/interface/EcalNextToDeadChannelRcd.h"
00008 #include "FWCore/Framework/interface/ESHandle.h"
00009 #include "FWCore/Framework/interface/EventSetup.h"
00010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00011 
00012 float EcalTools::swissCross( const DetId& id, 
00013                              const EcalRecHitCollection & recHits, 
00014                              float recHitThreshold , 
00015                              bool avoidIeta85){
00016   // compute swissCross
00017   if ( id.subdetId() == EcalBarrel ) {
00018     EBDetId ebId( id );
00019     // avoid recHits at |eta|=85 where one side of the neighbours is missing
00020     // (may improve considering also eta module borders, but no
00021     // evidence for the time being that there the performance is
00022     // different)
00023     if ( abs(ebId.ieta())==85 && avoidIeta85) return 0;
00024     // select recHits with Et above recHitThreshold
00025     if ( recHitApproxEt( id, recHits ) < recHitThreshold ) return 0;
00026     float s4 = 0;
00027     float e1 = recHitE( id, recHits );
00028     // protect against nan (if 0 threshold is given above)
00029     if ( e1 == 0 ) return 0;
00030     s4 += recHitE( id, recHits,  1,  0 );
00031     s4 += recHitE( id, recHits, -1,  0 );
00032     s4 += recHitE( id, recHits,  0,  1 );
00033     s4 += recHitE( id, recHits,  0, -1 );
00034     return 1 - s4 / e1;
00035   } else if ( id.subdetId() == EcalEndcap ) {
00036     EEDetId eeId( id );
00037     // select recHits with E above recHitThreshold
00038     float e1 = recHitE( id, recHits );
00039     if ( e1 < recHitThreshold ) return 0;
00040     float s4 = 0;
00041     // protect against nan (if 0 threshold is given above)
00042     if ( e1 == 0 ) return 0;
00043     s4 += recHitE( id, recHits,  1,  0 );
00044     s4 += recHitE( id, recHits, -1,  0 );
00045     s4 += recHitE( id, recHits,  0,  1 );
00046     s4 += recHitE( id, recHits,  0, -1 );
00047     return 1 - s4 / e1;
00048   }
00049   return 0;
00050 }
00051 
00052 
00053 bool EcalTools::isNextToDead( const DetId& id, const edm::EventSetup& es){
00054   
00055   edm::ESHandle<EcalNextToDeadChannel> dch;
00056   es.get<EcalNextToDeadChannelRcd>().get(dch);
00057   
00058   EcalNextToDeadChannel::const_iterator chIt = dch->find( id );
00059 
00060   if ( chIt != dch->end() ) {
00061     return *chIt;
00062   } else {
00063     edm::LogError("EcalDBError") 
00064       << "No NextToDead  status found for xtal " 
00065       << id.rawId() ;
00066   }
00067 
00068   return false;
00069   
00070 }
00071 
00072 bool EcalTools::isNextToDeadFromNeighbours( const DetId& id, 
00073                                             const EcalChannelStatus& chs,
00074                                             int chStatusThreshold) {
00075   
00076   if (deadNeighbour(id,chs, chStatusThreshold, 1, 0)) return true;
00077   if (deadNeighbour(id,chs, chStatusThreshold,-1, 0)) return true;
00078   if (deadNeighbour(id,chs, chStatusThreshold, 0, 1)) return true;
00079   if (deadNeighbour(id,chs, chStatusThreshold, 0,-1)) return true;
00080   if (deadNeighbour(id,chs, chStatusThreshold, 1,-1)) return true;
00081   if (deadNeighbour(id,chs, chStatusThreshold, 1, 1)) return true;
00082   if (deadNeighbour(id,chs, chStatusThreshold,-1, 1)) return true;
00083   if (deadNeighbour(id,chs, chStatusThreshold,-1,-1)) return true;
00084 
00085   return false;
00086 }
00087 
00088 
00089 bool EcalTools::deadNeighbour(const DetId& id, 
00090                               const EcalChannelStatus& chs,
00091                               int chStatusThreshold, 
00092                               int dx, int dy){
00093 
00094   
00095   DetId nid;
00096   if( id.subdetId() == EcalBarrel) nid = EBDetId::offsetBy( id, dx, dy );
00097   else if( id.subdetId() == EcalEndcap) nid = EEDetId::offsetBy( id, dx, dy );
00098 
00099   if (!nid) return false;
00100 
00101   EcalChannelStatus::const_iterator chIt = chs.find( nid );
00102   uint16_t dbStatus = 0;
00103   if ( chIt != chs.end() ) {
00104     // note that 
00105     dbStatus = chIt->getStatusCode() ;
00106   } else {
00107     edm::LogError("EcalDBError") 
00108       << "No channel status found for xtal " 
00109       << id.rawId() 
00110       << "! something wrong with EcalChannelStatus in your DB? ";
00111   }
00112 
00113   return (dbStatus>=chStatusThreshold );
00114   
00115 }
00116 
00117 
00118 
00119 float EcalTools::recHitE( const DetId id, 
00120                           const EcalRecHitCollection & recHits,
00121                           int di, int dj )
00122 {
00123   // in the barrel:   di = dEta   dj = dPhi
00124   // in the endcap:   di = dX     dj = dY
00125   
00126   DetId nid;
00127   if( id.subdetId() == EcalBarrel) nid = EBDetId::offsetBy( id, di, dj );
00128   else if( id.subdetId() == EcalEndcap) nid = EEDetId::offsetBy( id, di, dj );
00129   
00130   return ( nid == DetId(0) ? 0 : recHitE( nid, recHits ) );
00131 }
00132 
00133 float EcalTools::recHitE( const DetId id, const EcalRecHitCollection &recHits ){
00134   if ( id == DetId(0) ) {
00135     return 0;
00136   } else {
00137     EcalRecHitCollection::const_iterator it = recHits.find( id );
00138     if ( it != recHits.end() ) return (*it).energy();
00139   }
00140   return 0;
00141 }
00142 
00143 float EcalTools::recHitApproxEt( const DetId id, const EcalRecHitCollection &recHits ){
00144   // for the time being works only for the barrel
00145   if ( id.subdetId() == EcalBarrel ) {
00146     return recHitE( id, recHits ) / cosh( EBDetId::approxEta( id ) );
00147   }
00148   return 0;
00149 }
00150 
00151 
00152 bool isNextToBoundary(const DetId& id){
00153 
00154   if ( id.subdetId() == EcalBarrel ) 
00155     return EBDetId::isNextToBoundary(EBDetId(id));
00156   else  if ( id.subdetId() == EcalEndcap )
00157     return EEDetId::isNextToBoundary(EEDetId(id));
00158 
00159   return false;
00160 }