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