CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_0/src/RecoEgamma/EgammaIsolationAlgos/plugins/EgammaRecHitExtractor.cc

Go to the documentation of this file.
00001 //*****************************************************************************
00002 // File:      EgammaRecHitExtractor.cc
00003 // ----------------------------------------------------------------------------
00004 // OrigAuth:  Matthias Mozer, hacked by Sam Harper (ie the ugly stuff is mine)
00005 // Institute: IIHE-VUB, RAL
00006 //=============================================================================
00007 //*****************************************************************************
00008 //C++ includes
00009 #include <vector>
00010 #include <functional>
00011 
00012 //ROOT includes
00013 #include <Math/VectorUtil.h>
00014 
00015 //CMSSW includes
00016 #include "CommonTools/Utils/interface/StringToEnumValue.h"
00017 
00018 #include "RecoEgamma/EgammaIsolationAlgos/plugins/EgammaRecHitExtractor.h"
00019 #include "DataFormats/DetId/interface/DetId.h"
00020 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00021 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00022 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00023 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00024 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00025 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00026 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
00027 #include "RecoCaloTools/MetaCollections/interface/CaloRecHitMetaCollections.h"
00028 #include "DataFormats/Math/interface/deltaPhi.h"
00029 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
00030 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
00031 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00032 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
00033 
00034 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
00035 
00036 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00037 
00038 using namespace std;
00039 using namespace egammaisolation;
00040 using namespace reco::isodeposit;
00041 
00042 EgammaRecHitExtractor::EgammaRecHitExtractor(const edm::ParameterSet& par) : 
00043     etMin_(par.getParameter<double>("etMin")),
00044     energyMin_(par.getParameter<double>("energyMin")),
00045     extRadius_(par.getParameter<double>("extRadius")),
00046     intRadius_(par.getParameter<double>("intRadius")),
00047     intStrip_(par.getParameter<double>("intStrip")),
00048     barrelEcalHitsTag_(par.getParameter<edm::InputTag>("barrelEcalHits")), 
00049     endcapEcalHitsTag_(par.getParameter<edm::InputTag>("endcapEcalHits")),
00050     fakeNegativeDeposit_(par.getParameter<bool>("subtractSuperClusterEnergy")),
00051     tryBoth_(par.getParameter<bool>("tryBoth")),
00052     vetoClustered_(par.getParameter<bool>("vetoClustered")),
00053     sameTag_(false),
00054     severityLevelCut_(par.getParameter<int>("severityLevelCut"))
00055     //severityRecHitThreshold_(par.getParameter<double>("severityRecHitThreshold")),
00056     //spIdString_(par.getParameter<std::string>("spikeIdString")),
00057     //spIdThreshold_(par.getParameter<double>("spikeIdThreshold")),
00058 
00059 { 
00060   const std::vector<std::string> flagnames = 
00061     par.getParameter<std::vector<std::string> >("recHitFlagsToBeExcluded");
00062 
00063   v_chstatus_= StringToEnumValue<EcalRecHit::Flags>(flagnames);
00064 
00065 //     if     ( !spIdString_.compare("kE1OverE9") )                  spId_ = EcalSeverityLevelAlgo::kE1OverE9;
00066 //     else if( !spIdString_.compare("kSwissCross") )                spId_ = EcalSeverityLevelAlgo::kSwissCross;
00067 //     else if( !spIdString_.compare("kSwissCrossBordersIncluded") ) spId_ = EcalSeverityLevelAlgo::kSwissCrossBordersIncluded;
00068 //     else                                                          spId_ = EcalSeverityLevelAlgo::kSwissCross;
00069 
00070     if ((intRadius_ != 0.0) && (fakeNegativeDeposit_)) {
00071         throw cms::Exception("Configuration Error") << "EgammaRecHitExtractor: " << 
00072             "If you use 'subtractSuperClusterEnergy', you *must* set 'intRadius' to ZERO; it does not make sense, otherwise.";
00073     }
00074     std::string isoVariable = par.getParameter<std::string>("isolationVariable");
00075     if (isoVariable == "et") {
00076         useEt_ = true;
00077     } else if (isoVariable == "energy") {
00078         useEt_ = false;
00079     } else {
00080         throw cms::Exception("Configuration Error") << "EgammaRecHitExtractor: isolationVariable '" << isoVariable << "' not known. " 
00081             << " Supported values are 'et', 'energy'. ";
00082     }
00083     if (endcapEcalHitsTag_.encode() ==  barrelEcalHitsTag_.encode()) {
00084         sameTag_ = true;
00085         if (tryBoth_) {
00086             edm::LogWarning("EgammaRecHitExtractor") << "If you have configured 'barrelRecHits' == 'endcapRecHits', so I'm switching 'tryBoth' to FALSE.";
00087             tryBoth_ = false;
00088         }
00089     }
00090 
00091 }
00092 
00093 EgammaRecHitExtractor::~EgammaRecHitExtractor() { }
00094 
00095 reco::IsoDeposit EgammaRecHitExtractor::deposit(const edm::Event & iEvent, 
00096         const edm::EventSetup & iSetup, const reco::Candidate &emObject ) const {
00097     edm::ESHandle<CaloGeometry> pG;
00098     iSetup.get<CaloGeometryRecord>().get(pG);
00099 
00100     //Get the channel status from the db
00101     edm::ESHandle<EcalChannelStatus> chStatus;
00102     iSetup.get<EcalChannelStatusRcd>().get(chStatus);
00103 
00104     edm::ESHandle<EcalSeverityLevelAlgo> sevlv;
00105     iSetup.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
00106 
00107     const CaloGeometry* caloGeom = pG.product(); 
00108     const CaloSubdetectorGeometry* barrelgeom = caloGeom->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
00109     const CaloSubdetectorGeometry* endcapgeom = caloGeom->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
00110 
00111     static std::string metname = "EgammaIsolationAlgos|EgammaRecHitExtractor";
00112 
00113     std::auto_ptr<const CaloRecHitMetaCollectionV> barrelRecHits(0), endcapRecHits(0);
00114 
00115     //Get barrel ECAL RecHits 
00116     edm::Handle<EcalRecHitCollection> barrelEcalRecHitsH;
00117     iEvent.getByLabel(barrelEcalHitsTag_, barrelEcalRecHitsH);
00118 
00119     //Get endcap ECAL RecHits 
00120     edm::Handle<EcalRecHitCollection> endcapEcalRecHitsH;
00121     iEvent.getByLabel(endcapEcalHitsTag_, endcapEcalRecHitsH);
00122 
00123     //define isodeposit starting from candidate
00124     reco::SuperClusterRef sc = emObject.get<reco::SuperClusterRef>();
00125     math::XYZPoint caloPosition = sc->position();
00126 
00127     Direction candDir(caloPosition.eta(), caloPosition.phi());
00128     reco::IsoDeposit deposit( candDir );
00129     deposit.setVeto( reco::IsoDeposit::Veto(candDir, intRadius_) ); 
00130     double sinTheta = sin(2*atan(exp(-sc->eta())));
00131     deposit.addCandEnergy(sc->energy() * (useEt_ ? sinTheta : 1.0)) ;
00132 
00133     // subtract supercluster if desired
00134     double fakeEnergy = -sc->rawEnergy();
00135     if (fakeNegativeDeposit_) {
00136         deposit.addDeposit(candDir, fakeEnergy * (useEt_ ?  sinTheta : 1.0)); // not exactly clean...
00137     }
00138 
00139     // fill rechits
00140     bool inBarrel = sameTag_ || ( abs(sc->eta()) < 1.479 ); //check for barrel. If only one collection is used, use barrel
00141     if (inBarrel || tryBoth_) {
00142         collect(deposit, sc, barrelgeom, caloGeom, *barrelEcalRecHitsH, chStatus.product(), sevlv.product(), true);
00143     } 
00144     if ((!inBarrel) || tryBoth_) {
00145       collect(deposit, sc, endcapgeom, caloGeom, *endcapEcalRecHitsH, chStatus.product(), sevlv.product(), false);
00146     }
00147 
00148     return deposit;
00149 }
00150 
00151 void EgammaRecHitExtractor::collect(reco::IsoDeposit &deposit, 
00152                                     const reco::SuperClusterRef& sc, const CaloSubdetectorGeometry* subdet, 
00153                                     const CaloGeometry* caloGeom,
00154                                     const EcalRecHitCollection &hits, 
00155                                     const EcalChannelStatus* chStatus,
00156                                     const EcalSeverityLevelAlgo* sevLevel, 
00157                                     bool barrel) const 
00158 {
00159 
00160     GlobalPoint caloPosition(sc->position().x(), sc->position().y() , sc->position().z());
00161     CaloSubdetectorGeometry::DetIdSet chosen = subdet->getCells(caloPosition,extRadius_);
00162     EcalRecHitCollection::const_iterator j=hits.end();
00163     double caloeta=caloPosition.eta();
00164     double calophi=caloPosition.phi();
00165     double r2 = intRadius_*intRadius_;
00166 
00167     std::vector< std::pair<DetId, float> >::const_iterator rhIt;
00168 
00169 
00170     for (CaloSubdetectorGeometry::DetIdSet::const_iterator i = chosen.begin(), end = chosen.end() ; i != end;  ++i)  {  
00171         j=hits.find(*i);
00172         if(j != hits.end()){
00173             const  GlobalPoint & position = caloGeom->getPosition(*i);
00174             double eta = position.eta();
00175             double phi = position.phi();
00176             double energy = j->energy();
00177             double et = energy*position.perp()/position.mag();
00178             double phiDiff= reco::deltaPhi(phi,calophi);
00179 
00180             //check if we are supposed to veto clustered and then do so
00181             if(vetoClustered_) {
00182 
00183                 //Loop over basic clusters:
00184                 bool isClustered = false;
00185                 for(    reco::CaloCluster_iterator bcIt = sc->clustersBegin();bcIt != sc->clustersEnd(); ++bcIt) {
00186                     for(rhIt = (*bcIt)->hitsAndFractions().begin();rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
00187                         if( rhIt->first == *i ) isClustered = true;
00188                         if( isClustered ) break;
00189                     }
00190                     if( isClustered ) break;
00191                 } //end loop over basic clusters
00192 
00193                 if(isClustered) continue;
00194             }  //end if removeClustered
00195 
00196             //make sure we have a barrel rechit                                     
00197             //call the severity level method                                        
00198             //passing the EBDetId                                                   
00199             //the rechit collection in order to calculate the swiss crss            
00200             //and the EcalChannelRecHitRcd                                          
00201             //only consider rechits with ET >                                       
00202             //the SpikeId method (currently kE1OverE9 or kSwissCross)               
00203             //cut value for above                                                   
00204             //then if the severity level is too high, we continue to the next rechit
00205             if(barrel && sevLevel->severityLevel(EBDetId(j->id()), hits) >= severityLevelCut_)
00206               continue;                                 
00207             //                   *chStatus,                            
00208             //       severityRecHitThreshold_,             
00209             //       spId_,                                
00210             //       spIdThreshold_                        
00211             //   ) >= severityLevelCut_) continue;         
00212 
00213             //Check based on flags to protect from recovered channels from non-read towers
00214             //Assumption is that v_chstatus_ is empty unless doFlagChecks() has been called
00215             std::vector<int>::const_iterator vit = std::find( v_chstatus_.begin(), v_chstatus_.end(),  ((EcalRecHit*)(&*j))->recoFlag() );
00216             if ( vit != v_chstatus_.end() ) continue; // the recHit has to be excluded from the iso sum
00217 
00218 
00219             if(et > etMin_ 
00220                && energy > energyMin_  //Changed to fabs - then changed back to energy
00221                && fabs(eta-caloeta) > intStrip_ 
00222                && (eta-caloeta)*(eta-caloeta) + phiDiff*phiDiff >r2 ) {
00223               
00224                 deposit.addDeposit( Direction(eta, phi), (useEt_ ? et : energy) );
00225 
00226             }
00227         }
00228     }
00229 } 
00230 
00231