00001 #include "RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.h"
00002
00003 #include "FWCore/Framework/interface/EventSetup.h"
00004 #include "FWCore/Framework/interface/Event.h"
00005 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
00006 #include "CondFormats/DataRecord/interface/EcalTimeCalibConstantsRcd.h"
00007 #include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"
00008 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
00009 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
00010 #include "DataFormats/EcalDetId/interface/EcalScDetId.h"
00011 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
00012 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00013 #include "Geometry/EcalMapping/interface/EcalMappingRcd.h"
00014
00015 #include "CondFormats/EcalObjects/interface/EcalTimeCalibConstants.h"
00016 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbRecord.h"
00017
00018 #include "RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryAlgos.h"
00019
00020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00021
00022 EcalRecHitWorkerRecover::EcalRecHitWorkerRecover(const edm::ParameterSet&ps) :
00023 EcalRecHitWorkerBaseClass(ps)
00024 {
00025 rechitMaker_ = new EcalRecHitSimpleAlgo();
00026
00027 singleRecoveryMethod_ = ps.getParameter<std::string>("singleChannelRecoveryMethod");
00028 singleRecoveryThreshold_ = ps.getParameter<double>("singleChannelRecoveryThreshold");
00029 killDeadChannels_ = ps.getParameter<bool>("killDeadChannels");
00030 recoverEBIsolatedChannels_ = ps.getParameter<bool>("recoverEBIsolatedChannels");
00031 recoverEEIsolatedChannels_ = ps.getParameter<bool>("recoverEEIsolatedChannels");
00032 recoverEBVFE_ = ps.getParameter<bool>("recoverEBVFE");
00033 recoverEEVFE_ = ps.getParameter<bool>("recoverEEVFE");
00034 recoverEBFE_ = ps.getParameter<bool>("recoverEBFE");
00035 recoverEEFE_ = ps.getParameter<bool>("recoverEEFE");
00036
00037 tpDigiCollection_ = ps.getParameter<edm::InputTag>("triggerPrimitiveDigiCollection");
00038 logWarningEtThreshold_EB_FE_ = ps.getParameter<double>("logWarningEtThreshold_EB_FE");
00039 logWarningEtThreshold_EE_FE_ = ps.getParameter<double>("logWarningEtThreshold_EE_FE");
00040 }
00041
00042
00043 void EcalRecHitWorkerRecover::set(const edm::EventSetup& es)
00044 {
00045
00046 es.get<EcalLaserDbRecord>().get(laser);
00047 es.get<CaloTopologyRecord>().get(caloTopology_);
00048 ecalScale_.setEventSetup( es );
00049 es.get<EcalMappingRcd>().get(pEcalMapping_);
00050 ecalMapping_ = pEcalMapping_.product();
00051
00052 es.get<EcalBarrelGeometryRecord>().get("EcalBarrel",pEBGeom_);
00053 es.get<EcalEndcapGeometryRecord>().get("EcalEndcap",pEEGeom_);
00054 es.get<CaloGeometryRecord>().get(caloGeometry_);
00055 geo_ = caloGeometry_.product();
00056 ebGeom_ = pEBGeom_.product();
00057 eeGeom_ = pEEGeom_.product();
00058 es.get<IdealGeometryRecord>().get(ttMap_);
00059 recoveredDetIds_EB_.clear();
00060 recoveredDetIds_EE_.clear();
00061 }
00062
00063
00064 bool
00065 EcalRecHitWorkerRecover::run( const edm::Event & evt,
00066 const EcalUncalibratedRecHit& uncalibRH,
00067 EcalRecHitCollection & result )
00068 {
00069 DetId detId=uncalibRH.id();
00070 uint32_t flags = uncalibRH.recoFlag();
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081 if ( killDeadChannels_ ) {
00082 if ( (flags == EcalRecHitWorkerRecover::EB_single && !recoverEBIsolatedChannels_)
00083 || (flags == EcalRecHitWorkerRecover::EE_single && !recoverEEIsolatedChannels_)
00084 || (flags == EcalRecHitWorkerRecover::EB_VFE && !recoverEBVFE_)
00085 || (flags == EcalRecHitWorkerRecover::EE_VFE && !recoverEEVFE_)
00086 ) {
00087 EcalRecHit hit( detId, 0., 0., EcalRecHit::kDead );
00088 hit.setFlagBits( (0x1 << EcalRecHit::kDead) ) ;
00089 insertRecHit( hit, result);
00090 return true;
00091 }
00092 if ( flags == EcalRecHitWorkerRecover::EB_FE && !recoverEBFE_) {
00093 EcalTrigTowerDetId ttDetId( ((EBDetId)detId).tower() );
00094 std::vector<DetId> vid = ttMap_->constituentsOf( ttDetId );
00095 for ( std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); ++dit ) {
00096 EcalRecHit hit( (*dit), 0., 0., EcalRecHit::kDead );
00097 hit.setFlagBits( (0x1 << EcalRecHit::kDead) ) ;
00098 insertRecHit( hit, result );
00099 }
00100 if(logWarningEtThreshold_EB_FE_<0)return true;
00101 }
00102 if ( flags == EcalRecHitWorkerRecover::EE_FE && !recoverEEFE_) {
00103 EEDetId id( detId );
00104 EcalScDetId sc( 1+(id.ix()-1)/5, 1+(id.iy()-1)/5, id.zside() );
00105 std::vector<DetId> eeC;
00106 for(int dx=1; dx<=5; ++dx){
00107 for(int dy=1; dy<=5; ++dy){
00108 int ix = (sc.ix()-1)*5 + dx;
00109 int iy = (sc.iy()-1)*5 + dy;
00110 int iz = sc.zside();
00111 if(EEDetId::validDetId(ix, iy, iz)){
00112 eeC.push_back(EEDetId(ix, iy, iz));
00113 }
00114 }
00115 }
00116 for ( size_t i = 0; i < eeC.size(); ++i ) {
00117 EcalRecHit hit( eeC[i], 0., 0., EcalRecHit::kDead );
00118 hit.setFlagBits( (0x1 << EcalRecHit::kDead) ) ;
00119 insertRecHit( hit, result );
00120 }
00121 if(logWarningEtThreshold_EE_FE_<0) return true;
00122 }
00123 }
00124
00125 if ( flags == EcalRecHitWorkerRecover::EB_single ) {
00126
00127 const EcalRecHitCollection * hit_collection = &result;
00128 EcalDeadChannelRecoveryAlgos deadChannelCorrector(caloTopology_.product());
00129
00130
00131 EcalRecHit hit = deadChannelCorrector.correct( detId, hit_collection, singleRecoveryMethod_, singleRecoveryThreshold_ );
00132 EcalRecHitCollection::const_iterator ti = result.find( detId );
00133 if ( hit.energy() != 0 ) {
00134 hit.setFlags( EcalRecHit::kNeighboursRecovered );
00135 hit.setFlagBits( (0x1 << EcalRecHit::kNeighboursRecovered) ) ;
00136 } else {
00137
00138 hit.setFlags( EcalRecHit::kDead );
00139 hit.setFlagBits( (0x1 << EcalRecHit::kDead) ) ;
00140 }
00141 insertRecHit( hit, result );
00142 } else if ( flags == EcalRecHitWorkerRecover::EB_VFE ) {
00143
00144 EcalRecHit hit( detId, 0., 0., EcalRecHit::kDead );
00145 hit.setFlagBits( (0x1 << EcalRecHit::kDead) ) ;
00146
00147 insertRecHit( hit, result );
00148 } else if ( flags == EcalRecHitWorkerRecover::EB_FE ) {
00149
00150 EcalTrigTowerDetId ttDetId( ((EBDetId)detId).tower() );
00151 edm::Handle<EcalTrigPrimDigiCollection> pTPDigis;
00152 evt.getByLabel(tpDigiCollection_, pTPDigis);
00153 const EcalTrigPrimDigiCollection * tpDigis = 0;
00154 if ( pTPDigis.isValid() ) {
00155 tpDigis = pTPDigis.product();
00156 } else {
00157 edm::LogError("EcalRecHitWorkerRecover") << "Can't get the product " << tpDigiCollection_.instance()
00158 << " with label " << tpDigiCollection_.label();
00159 return false;
00160 }
00161 EcalTrigPrimDigiCollection::const_iterator tp = tpDigis->find( ttDetId );
00162
00163 if ( tp != tpDigis->end() ) {
00164
00165 std::vector<DetId> vid = ttMap_->constituentsOf( ttDetId );
00166 float tpEt = ecalScale_.getTPGInGeV( tp->compressedEt(), tp->id() );
00167 float tpEtThreshEB = logWarningEtThreshold_EB_FE_;
00168 if(tpEt>tpEtThreshEB){
00169 edm::LogWarning("EnergyInDeadEB_FE")<<"TP energy in the dead TT = "<<tpEt<<" at "<<ttDetId;
00170 }
00171 if ( !killDeadChannels_ || recoverEBFE_ ) {
00172
00173 for ( std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); ++dit ) {
00174 if (alreadyInserted(*dit)) continue;
00175 float theta = ebGeom_->getGeometry(*dit)->getPosition().theta();
00176 float tpEt = ecalScale_.getTPGInGeV( tp->compressedEt(), tp->id() );
00177 EcalRecHit hit( *dit, tpEt / (float)vid.size() / sin(theta), 0., EcalRecHit::kTowerRecovered );
00178 hit.setFlagBits( (0x1 << EcalRecHit::kTowerRecovered) ) ;
00179 if ( tp->compressedEt() == 0xFF ) hit.setFlagBits( (0x1 << EcalRecHit::kTPSaturated) );
00180 if ( tp->sFGVB() ) hit.setFlagBits( (0x1 << EcalRecHit::kL1SpikeFlag) );
00181 insertRecHit( hit, result );
00182 }
00183 }
00184 } else {
00185
00186 std::vector<DetId> vid = ttMap_->constituentsOf( ttDetId );
00187 for ( std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); ++dit ) {
00188 if (alreadyInserted(*dit)) continue;
00189 EcalRecHit hit( *dit,0., 0., EcalRecHit::kDead );
00190 hit.setFlagBits( (0x1 << EcalRecHit::kDead) ) ;
00191 EcalRecHitCollection::iterator it = result.find( *dit );
00192 insertRecHit( hit, result );
00193 }
00194 }
00195 } else if ( flags == EcalRecHitWorkerRecover::EE_FE ) {
00196
00197
00198
00199
00200
00201
00202
00203
00204 EEDetId eeId( detId );
00205 EcalScDetId sc( (eeId.ix()-1)/5+1, (eeId.iy()-1)/5+1, eeId.zside() );
00206 std::set<DetId> eeC;
00207 for(int dx=1; dx<=5; ++dx){
00208 for(int dy=1; dy<=5; ++dy){
00209 int ix = (sc.ix()-1)*5 + dx;
00210 int iy = (sc.iy()-1)*5 + dy;
00211 int iz = sc.zside();
00212 if(EEDetId::validDetId(ix, iy, iz)){
00213 eeC.insert(EEDetId(ix, iy, iz));
00214 }
00215 }
00216 }
00217
00218 edm::Handle<EcalTrigPrimDigiCollection> pTPDigis;
00219 evt.getByLabel(tpDigiCollection_, pTPDigis);
00220 const EcalTrigPrimDigiCollection * tpDigis = 0;
00221 if ( pTPDigis.isValid() ) {
00222 tpDigis = pTPDigis.product();
00223 } else {
00224 edm::LogError("EcalRecHitWorkerRecover") << "Can't get the product " << tpDigiCollection_.instance()
00225 << " with label " << tpDigiCollection_.label();
00226 return false;
00227 }
00228
00229 std::set<EcalTrigTowerDetId> aTT;
00230 for ( std::set<DetId>::const_iterator it = eeC.begin(); it!=eeC.end(); ++it ) {
00231 aTT.insert( ttMap_->towerOf( *it ) );
00232 }
00233
00234 float totE = 0;
00235
00236 std::set<DetId> aTTC;
00237 bool atLeastOneTPSaturated = false;
00238 for ( std::set<EcalTrigTowerDetId>::const_iterator it = aTT.begin(); it != aTT.end(); ++it ) {
00239
00240 EcalTrigPrimDigiCollection::const_iterator itTP = tpDigis->find( *it );
00241 if ( itTP != tpDigis->end() ) {
00242 EcalTrigTowerDetId ttId = itTP->id();
00243
00244 std::vector<DetId> v = ttMap_->constituentsOf( *it );
00245 if ( itTP->compressedEt() == 0xFF ){
00246 atLeastOneTPSaturated = true;
00247
00248
00249
00250
00251
00252 totE += estimateEnergy(itTP->id().ietaAbs(), &result, eeC, v);
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266 }
00267 else {totE += ((it->ietaAbs()>26)?2:1)*ecalScale_.getTPGInGeV( itTP->compressedEt(), itTP->id() );}
00268
00269
00270
00271
00272 if (itTP->compressedEt() == 0){
00273 for (size_t i = 0 ; i < v.size(); ++i)
00274 eeC.erase(v[i]);
00275 }
00276 else if (itTP->compressedEt()!=0xFF){
00277 for ( size_t j = 0; j < v.size(); ++j ) {
00278 aTTC.insert( v[j] );
00279 }
00280 }
00281
00282 }
00283 }
00284
00285
00286
00287
00288 for ( std::set<DetId>::const_iterator it = eeC.begin(); it != eeC.end(); ++it ) {
00289 aTTC.erase(*it);
00290 }
00291
00292 const EcalRecHitCollection * hits = &result;
00293 for ( std::set<DetId>::const_iterator it = aTTC.begin(); it != aTTC.end(); ++it ) {
00294 EcalRecHitCollection::const_iterator jt = hits->find( *it );
00295 if ( jt != hits->end() ) {
00296 float energy = jt->energy();
00297 float eta = geo_->getPosition(jt->id()).eta();
00298 float pf = 1.0/cosh(eta);
00299
00300
00301 totE -= energy*pf;
00302 }
00303 }
00304
00305
00306 float scEt = totE;
00307 float scEtThreshEE = logWarningEtThreshold_EE_FE_;
00308 if(scEt>scEtThreshEE){
00309 edm::LogWarning("EnergyInDeadEE_FE")<<"TP energy in the dead TT = "<<scEt<<" at "<<sc;
00310 }
00311
00312
00313 if ( !killDeadChannels_ || recoverEEFE_ ) {
00314
00315 for ( std::set<DetId>::const_iterator it = eeC.begin(); it != eeC.end(); ++it ) {
00316 EcalRecHit hit( *it, 0., 0., EcalRecHit::kDead );
00317 hit.setFlagBits( (0x1 << EcalRecHit::kDead) ) ;
00318 float eta = geo_->getPosition(*it).eta();
00319 float pf = 1.0/cosh(eta);
00320 hit = EcalRecHit( *it, totE / ((float)eeC.size()*pf), 0., EcalRecHit::kTowerRecovered );
00321 if (atLeastOneTPSaturated) hit.setFlagBits( (0x1 << EcalRecHit::kTPSaturated) );
00322
00323 insertRecHit( hit, result );
00324 }
00325 }
00326 }
00327 return true;
00328 }
00329
00330 float EcalRecHitWorkerRecover::estimateEnergy(int ieta, EcalRecHitCollection* hits, std::set<DetId> sId, std::vector<DetId> vId ) {
00331
00332 float xtalE=0;
00333 int count = 0;
00334 for (std::vector<DetId>::const_iterator vIdit = vId.begin(); vIdit != vId.end(); ++ vIdit){
00335 std::set<DetId>::const_iterator sIdit = sId.find(*vIdit);
00336 if (sIdit==sId.end()){
00337 float energy = hits->find(*vIdit)->energy();
00338 float eta = geo_->getPosition(*vIdit).eta();
00339 float pf = 1.0/cosh(eta);
00340 xtalE += energy*pf;
00341 count++;
00342 }
00343 }
00344
00345 if (count==0) return 63.75*(ieta>26?2:1);
00346 else return xtalE*((vId.size()/(float)count) - 1)*(ieta>26?2:1);
00347
00348
00349 }
00350
00351
00352 void EcalRecHitWorkerRecover::insertRecHit( const EcalRecHit &hit, EcalRecHitCollection &collection )
00353 {
00354
00355 if ( alreadyInserted( hit.id() ) ) {
00356 edm::LogWarning("EcalRecHitWorkerRecover") << "DetId already recovered! Skipping...";
00357 return;
00358 }
00359 EcalRecHitCollection::iterator it = collection.find( hit.id() );
00360 if ( it == collection.end() ) {
00361
00362 collection.push_back( hit );
00363 } else {
00364
00365 *it = hit;
00366 }
00367 if ( hit.id().subdetId() == EcalBarrel ) {
00368 recoveredDetIds_EB_.insert( hit.id() );
00369 } else if ( hit.id().subdetId() == EcalEndcap ) {
00370 recoveredDetIds_EE_.insert( hit.id() );
00371 } else {
00372 edm::LogError("EcalRecHitWorkerRecover::InvalidDetId") << "Invalid DetId " << hit.id().rawId();
00373 }
00374 }
00375
00376
00377 bool EcalRecHitWorkerRecover::alreadyInserted( const DetId & id )
00378 {
00379 bool res = false;
00380 if ( id.subdetId() == EcalBarrel ) {
00381 res = ( recoveredDetIds_EB_.find( id ) != recoveredDetIds_EB_.end() );
00382 } else if ( id.subdetId() == EcalEndcap ) {
00383 res = ( recoveredDetIds_EE_.find( id ) != recoveredDetIds_EE_.end() );
00384 } else {
00385 edm::LogError("EcalRecHitWorkerRecover::InvalidDetId") << "Invalid DetId " << id.rawId();
00386 }
00387 return res;
00388 }
00389
00390
00391
00392 float EcalRecHitWorkerRecover::recCheckCalib(float eTT, int ieta){
00393
00394 return eTT;
00395
00396 }
00397
00398
00399 #include "FWCore/Framework/interface/MakerMacros.h"
00400 #include "RecoLocalCalo/EcalRecProducers/interface/EcalRecHitWorkerFactory.h"
00401 DEFINE_EDM_PLUGIN( EcalRecHitWorkerFactory, EcalRecHitWorkerRecover, "EcalRecHitWorkerRecover" );