CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.cc

Go to the documentation of this file.
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         // isolated channel recovery
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         // geometry...
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         // get laser coefficient
00073         //float lasercalib = laser->getLaserCorrection( detId, evt.time());
00074 
00075         // killDeadChannels_ = true, means explicitely kill dead channels even if the recovered energies are computed in the code
00076         // if you don't want to store the recovered energies in the rechit you can produce LogWarnings if logWarningEtThreshold_EB(EE)_FE>0 
00077         // logWarningEtThreshold_EB(EE)_FE_<0 will not compute the recovered energies at all (faster)
00078         // Revovery in the EE is not tested, recovered energies may not make sense
00079         // EE recovery computation is not tested against segmentation faults, use with caution even if you are going to killDeadChannels=true
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); // insert trivial rechit with kDead flag
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 ); // insert trivial rechit with kDead flag
00099                         }
00100                         if(logWarningEtThreshold_EB_FE_<0)return true; // if you don't want log warning just 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 ); // insert trivial rechit with kDead flag
00120                         }
00121                         if(logWarningEtThreshold_EE_FE_<0)   return true; // if you don't want log warning just return true
00122                 }
00123         }
00124 
00125         if ( flags == EcalRecHitWorkerRecover::EB_single ) {
00126                 // recover as single dead channel
00127                 const EcalRecHitCollection * hit_collection = &result;
00128                 EcalDeadChannelRecoveryAlgos deadChannelCorrector(caloTopology_.product());
00129 
00130                 // channel recovery
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                         // recovery failed
00138                         hit.setFlags( EcalRecHit::kDead );
00139                         hit.setFlagBits( (0x1 << EcalRecHit::kDead) ) ;
00140                 }
00141                 insertRecHit( hit, result );
00142         } else if ( flags == EcalRecHitWorkerRecover::EB_VFE ) {
00143                 // recover as dead VFE
00144                 EcalRecHit hit( detId, 0., 0., EcalRecHit::kDead );
00145                 hit.setFlagBits( (0x1 << EcalRecHit::kDead) ) ;
00146                 // recovery not implemented
00147                 insertRecHit( hit, result );
00148         } else if ( flags == EcalRecHitWorkerRecover::EB_FE ) {
00149                 // recover as dead TT
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                 // recover the whole trigger tower
00163                 if ( tp != tpDigis->end() ) {
00164                         //std::vector<DetId> vid = ecalMapping_->dccTowerConstituents( ecalMapping_->DCCid( ttDetId ), ecalMapping_->iTT( ttDetId ) );
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                                 // democratic energy sharing
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                         // tp not found => recovery failed
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                         // Structure for recovery:
00197                         // ** SC --> EEDetId constituents (eeC) --> associated Trigger Towers (aTT) --> EEDetId constituents (aTTC)
00198                         // ** energy for a SC EEDetId = [ sum_aTT(energy) - sum_aTTC(energy) ] / N_eeC
00199                         // .. i.e. the total energy of the TTs covering the SC minus 
00200                         // .. the energy of the recHits in the TTs but not in the SC
00201                         //std::vector<DetId> vid = ecalMapping_->dccTowerConstituents( ecalMapping_->DCCid( ttDetId ), ecalMapping_->iTT( ttDetId ) );
00202                         // due to lack of implementation of the EcalTrigTowerDetId ix,iy methods in EE we compute Et recovered energies (in EB we compute E)
00203                         // --- RECOVERY NOT YET VALIDATED
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                         // associated trigger towers
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                         // associated trigger towers: total energy
00234                         float totE = 0;
00235                         // associated trigger towers: EEDetId constituents
00236                         std::set<DetId> aTTC;
00237                         bool atLeastOneTPSaturated = false;
00238                         for ( std::set<EcalTrigTowerDetId>::const_iterator it = aTT.begin(); it != aTT.end(); ++it ) {
00239                                 // add the energy of this trigger tower
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 ){ // In the case of a saturated trigger tower, a fraction
00246                                           atLeastOneTPSaturated = true; //of the saturated energy is put in: number of xtals in dead region/total xtals in TT *63.75
00247                                                 
00248                                           //Alternative recovery algorithm that I will now investigate.
00249                                           //Estimate energy sums the energy in the working channels, then decides how much energy
00250                                           //to put here depending on that. Duncan 20101203
00251                                           
00252                                           totE += estimateEnergy(itTP->id().ietaAbs(), &result, eeC, v);
00253                                           
00254                                           /* 
00255                                              These commented out lines use
00256                                              64GeV*fraction of the TT overlapping the dead FE
00257                                             
00258                                           int count = 0;
00259                                           for (std::vector<DetId>::const_iterator idsit = v.begin(); idsit != v.end(); ++ idsit){
00260                                           std::set<DetId>::const_iterator itFind = eeC.find(*idsit);
00261                                           if (itFind != eeC.end())
00262                                           ++count;
00263                                           }
00264                                           //std::cout << count << ", " << v.size() << std::endl;
00265                                           totE+=((float)count/(float)v.size())* ((it->ietaAbs()>26)?2*ecalScale_.getTPGInGeV( itTP->compressedEt(), itTP->id() ):ecalScale_.getTPGInGeV( itTP->compressedEt(), itTP->id() ));*/
00266                                         }
00267                                         else {totE += ((it->ietaAbs()>26)?2:1)*ecalScale_.getTPGInGeV( itTP->compressedEt(), itTP->id() );}
00268                                         
00269                                 
00270                                         // get the trigger tower constituents
00271                                         
00272                                         if (itTP->compressedEt() == 0){ // If there's no energy in TT, the constituents are removed from the recovery.
00273                                                  for (size_t i = 0 ; i < v.size(); ++i)
00274                                                          eeC.erase(v[i]);
00275                                         }
00276                                         else if (itTP->compressedEt()!=0xFF){ //If it's saturated the energy has already been determined, so we do not want to subtract any channels
00277                                          for ( size_t j = 0; j < v.size(); ++j ) {
00278                                                  aTTC.insert( v[j] );
00279                                          }
00280                                         }
00281                                     
00282                                 }
00283                         }
00284                         // remove crystals of dead SC
00285                         // (this step is not needed if sure that SC crystals are not 
00286                         // in the recHit collection)
00287                         
00288                         for ( std::set<DetId>::const_iterator it = eeC.begin(); it != eeC.end(); ++it ) {
00289                           aTTC.erase(*it);
00290                         }
00291                         // compute the total energy for the dead SC
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(); // Correct conversion to Et
00297                                   float eta = geo_->getPosition(jt->id()).eta();
00298                                   float pf = 1.0/cosh(eta);
00299                                   //float theta = eeGeom_->getGeometry( *it )->getPosition().theta();
00300                                   // use Et instead of E, consistent with the Et estimation of the associated TT
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                         // assign the energy to the SC crystals
00313                         if ( !killDeadChannels_ || recoverEEFE_ ) { // if eeC is empty, i.e. there are no hits 
00314                                                                     // in the tower, nothing is returned. No negative values from noise.
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(); //Convert back to E from Et for the recovered hits
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); //If there are no overlapping crystals return saturated value.
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         // skip already inserted DetId's and raise a log warning
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                 // insert the hit in the collection
00362                 collection.push_back( hit );
00363         } else {
00364                 // overwrite existing recHit
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 // In the future, this will be used to calibrate the TT energy. There is a dependance on
00391 // eta at lower energies that can be corrected for here after more validation.
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" );