00001 #include <iostream>
00002 #include <vector>
00003 #include <memory>
00004
00005
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 #include "FWCore/Utilities/interface/Exception.h"
00008
00009 #include "CommonTools/Utils/interface/StringToEnumValue.h"
00010
00011 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00012 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00013 #include "Geometry/CaloTopology/interface/CaloTopology.h"
00014 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
00015
00016
00017 #include "DataFormats/VertexReco/interface/Vertex.h"
00018 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00019 #include "DataFormats/EgammaReco/interface/ClusterShape.h"
00020 #include "DataFormats/EgammaCandidates/interface/PhotonCore.h"
00021 #include "DataFormats/EgammaCandidates/interface/Photon.h"
00022 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
00023 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
00024
00025 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
00026 #include "RecoCaloTools/Selectors/interface/CaloConeSelector.h"
00027
00028 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
00029 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00030
00031 #include "RecoEgamma/EgammaPhotonProducers/interface/GEDPhotonProducer.h"
00032 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTowerIsolation.h"
00033 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionBaseClass.h"
00034 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionFactory.h"
00035 #include "RecoEcal/EgammaCoreTools/plugins/EcalClusterCrackCorrection.h"
00036 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHadTower.h"
00037
00038 GEDPhotonProducer::GEDPhotonProducer(const edm::ParameterSet& config) :
00039
00040 conf_(config)
00041 {
00042
00043
00044
00045 photonCoreProducer_ = conf_.getParameter<edm::InputTag>("photonCoreProducer");
00046 barrelEcalHits_ = conf_.getParameter<edm::InputTag>("barrelEcalHits");
00047 endcapEcalHits_ = conf_.getParameter<edm::InputTag>("endcapEcalHits");
00048 vertexProducer_ = conf_.getParameter<std::string>("primaryVertexProducer");
00049 hcalTowers_ = conf_.getParameter<edm::InputTag>("hcalTowers");
00050 hOverEConeSize_ = conf_.getParameter<double>("hOverEConeSize");
00051 highEt_ = conf_.getParameter<double>("highEt");
00052
00053 minR9Barrel_ = conf_.getParameter<double>("minR9Barrel");
00054 minR9Endcap_ = conf_.getParameter<double>("minR9Endcap");
00055 usePrimaryVertex_ = conf_.getParameter<bool>("usePrimaryVertex");
00056 runMIPTagger_ = conf_.getParameter<bool>("runMIPTagger");
00057
00058 candidateP4type_ = config.getParameter<std::string>("candidateP4type") ;
00059
00060 edm::ParameterSet posCalcParameters =
00061 config.getParameter<edm::ParameterSet>("posCalcParameters");
00062 posCalculator_ = PositionCalc(posCalcParameters);
00063
00064
00065
00066
00067 const std::vector<std::string> flagnamesEB =
00068 config.getParameter<std::vector<std::string> >("RecHitFlagToBeExcludedEB");
00069
00070 const std::vector<std::string> flagnamesEE =
00071 config.getParameter<std::vector<std::string> >("RecHitFlagToBeExcludedEE");
00072
00073 flagsexclEB_=
00074 StringToEnumValue<EcalRecHit::Flags>(flagnamesEB);
00075
00076 flagsexclEE_=
00077 StringToEnumValue<EcalRecHit::Flags>(flagnamesEE);
00078
00079 const std::vector<std::string> severitynamesEB =
00080 config.getParameter<std::vector<std::string> >("RecHitSeverityToBeExcludedEB");
00081
00082 severitiesexclEB_=
00083 StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynamesEB);
00084
00085 const std::vector<std::string> severitynamesEE =
00086 config.getParameter<std::vector<std::string> >("RecHitSeverityToBeExcludedEE");
00087
00088 severitiesexclEE_=
00089 StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynamesEE);
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 preselCutValuesBarrel_.push_back(conf_.getParameter<double>("minSCEtBarrel"));
00106 preselCutValuesBarrel_.push_back(conf_.getParameter<double>("maxHoverEBarrel"));
00107 preselCutValuesBarrel_.push_back(conf_.getParameter<double>("ecalRecHitSumEtOffsetBarrel"));
00108 preselCutValuesBarrel_.push_back(conf_.getParameter<double>("ecalRecHitSumEtSlopeBarrel"));
00109 preselCutValuesBarrel_.push_back(conf_.getParameter<double>("hcalTowerSumEtOffsetBarrel"));
00110 preselCutValuesBarrel_.push_back(conf_.getParameter<double>("hcalTowerSumEtSlopeBarrel"));
00111 preselCutValuesBarrel_.push_back(conf_.getParameter<double>("nTrackSolidConeBarrel"));
00112 preselCutValuesBarrel_.push_back(conf_.getParameter<double>("nTrackHollowConeBarrel"));
00113 preselCutValuesBarrel_.push_back(conf_.getParameter<double>("trackPtSumSolidConeBarrel"));
00114 preselCutValuesBarrel_.push_back(conf_.getParameter<double>("trackPtSumHollowConeBarrel"));
00115 preselCutValuesBarrel_.push_back(conf_.getParameter<double>("sigmaIetaIetaCutBarrel"));
00116
00117 preselCutValuesEndcap_.push_back(conf_.getParameter<double>("minSCEtEndcap"));
00118 preselCutValuesEndcap_.push_back(conf_.getParameter<double>("maxHoverEEndcap"));
00119 preselCutValuesEndcap_.push_back(conf_.getParameter<double>("ecalRecHitSumEtOffsetEndcap"));
00120 preselCutValuesEndcap_.push_back(conf_.getParameter<double>("ecalRecHitSumEtSlopeEndcap"));
00121 preselCutValuesEndcap_.push_back(conf_.getParameter<double>("hcalTowerSumEtOffsetEndcap"));
00122 preselCutValuesEndcap_.push_back(conf_.getParameter<double>("hcalTowerSumEtSlopeEndcap"));
00123 preselCutValuesEndcap_.push_back(conf_.getParameter<double>("nTrackSolidConeEndcap"));
00124 preselCutValuesEndcap_.push_back(conf_.getParameter<double>("nTrackHollowConeEndcap"));
00125 preselCutValuesEndcap_.push_back(conf_.getParameter<double>("trackPtSumSolidConeEndcap"));
00126 preselCutValuesEndcap_.push_back(conf_.getParameter<double>("trackPtSumHollowConeEndcap"));
00127 preselCutValuesEndcap_.push_back(conf_.getParameter<double>("sigmaIetaIetaCutEndcap"));
00128
00129
00130
00131 produces< reco::PhotonCollection >(PhotonCollection_);
00132
00133 }
00134
00135 GEDPhotonProducer::~GEDPhotonProducer()
00136 {
00137
00138
00139 }
00140
00141
00142
00143 void GEDPhotonProducer::beginRun (edm::Run const& r, edm::EventSetup const & theEventSetup) {
00144
00145 thePhotonIsolationCalculator_ = new PhotonIsolationCalculator();
00146 edm::ParameterSet isolationSumsCalculatorSet = conf_.getParameter<edm::ParameterSet>("isolationSumsCalculatorSet");
00147 thePhotonIsolationCalculator_->setup(isolationSumsCalculatorSet, flagsexclEB_, flagsexclEE_, severitiesexclEB_, severitiesexclEE_);
00148
00149 thePhotonMIPHaloTagger_ = new PhotonMIPHaloTagger();
00150 edm::ParameterSet mipVariableSet = conf_.getParameter<edm::ParameterSet>("mipVariableSet");
00151 thePhotonMIPHaloTagger_->setup(mipVariableSet);
00152 thePhotonEnergyCorrector_ = new PhotonEnergyCorrector(conf_);
00153 thePhotonEnergyCorrector_ -> init(theEventSetup);
00154 }
00155
00156 void GEDPhotonProducer::endRun (edm::Run const& r, edm::EventSetup const & theEventSetup) {
00157
00158 delete thePhotonIsolationCalculator_;
00159 delete thePhotonMIPHaloTagger_;
00160 delete thePhotonEnergyCorrector_;
00161 }
00162
00163
00164 void GEDPhotonProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
00165
00166 using namespace edm;
00167
00168
00169 reco::PhotonCollection outputPhotonCollection;
00170 std::auto_ptr< reco::PhotonCollection > outputPhotonCollection_p(new reco::PhotonCollection);
00171
00172
00173
00174 bool validPhotonCoreHandle=true;
00175 Handle<reco::PhotonCoreCollection> photonCoreHandle;
00176 theEvent.getByLabel(photonCoreProducer_,photonCoreHandle);
00177 if (!photonCoreHandle.isValid()) {
00178 edm::LogError("GEDPhotonProducer") << "Error! Can't get the product "<<photonCoreProducer_.label();
00179 validPhotonCoreHandle=false;
00180 }
00181
00182
00183 bool validEcalRecHits=true;
00184 Handle<EcalRecHitCollection> barrelHitHandle;
00185 EcalRecHitCollection barrelRecHits;
00186 theEvent.getByLabel(barrelEcalHits_, barrelHitHandle);
00187 if (!barrelHitHandle.isValid()) {
00188 edm::LogError("GEDPhotonProducer") << "Error! Can't get the product "<<barrelEcalHits_.label();
00189 validEcalRecHits=false;
00190 }
00191 if ( validEcalRecHits) barrelRecHits = *(barrelHitHandle.product());
00192
00193
00194 Handle<EcalRecHitCollection> endcapHitHandle;
00195 theEvent.getByLabel(endcapEcalHits_, endcapHitHandle);
00196 EcalRecHitCollection endcapRecHits;
00197 if (!endcapHitHandle.isValid()) {
00198 edm::LogError("GEDPhotonProducer") << "Error! Can't get the product "<<endcapEcalHits_.label();
00199 validEcalRecHits=false;
00200 }
00201 if( validEcalRecHits) endcapRecHits = *(endcapHitHandle.product());
00202
00203
00204
00205 edm::ESHandle<EcalSeverityLevelAlgo> sevLv;
00206 theEventSetup.get<EcalSeverityLevelAlgoRcd>().get(sevLv);
00207
00208
00209
00210
00211 Handle<CaloTowerCollection> hcalTowersHandle;
00212 theEvent.getByLabel(hcalTowers_, hcalTowersHandle);
00213
00214
00215
00216 theEventSetup.get<CaloGeometryRecord>().get(theCaloGeom_);
00217
00218
00219
00220
00221
00222 edm::ESHandle<CaloTopology> pTopology;
00223 theEventSetup.get<CaloTopologyRecord>().get(theCaloTopo_);
00224 const CaloTopology *topology = theCaloTopo_.product();
00225
00226
00227 Handle<reco::VertexCollection> vertexHandle;
00228 reco::VertexCollection vertexCollection;
00229 bool validVertex=true;
00230 if ( usePrimaryVertex_ ) {
00231 theEvent.getByLabel(vertexProducer_, vertexHandle);
00232 if (!vertexHandle.isValid()) {
00233 edm::LogWarning("GEDPhotonProducer") << "Error! Can't get the product primary Vertex Collection "<< "\n";
00234 validVertex=false;
00235 }
00236 if (validVertex) vertexCollection = *(vertexHandle.product());
00237 }
00238
00239
00240
00241
00242 int iSC=0;
00243
00244 if ( validPhotonCoreHandle)
00245 fillPhotonCollection(theEvent,
00246 theEventSetup,
00247 photonCoreHandle,
00248 topology,
00249 &barrelRecHits,
00250 &endcapRecHits,
00251 hcalTowersHandle,
00252
00253 vertexCollection,
00254 outputPhotonCollection,
00255 iSC,
00256 sevLv.product());
00257
00258
00259
00260 edm::LogInfo("GEDPhotonProducer") << " Put in the event " << iSC << " Photon Candidates \n";
00261 outputPhotonCollection_p->assign(outputPhotonCollection.begin(),outputPhotonCollection.end());
00262 theEvent.put( outputPhotonCollection_p, PhotonCollection_);
00263
00264 }
00265
00266 void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt,
00267 edm::EventSetup const & es,
00268 const edm::Handle<reco::PhotonCoreCollection> & photonCoreHandle,
00269 const CaloTopology* topology,
00270 const EcalRecHitCollection* ecalBarrelHits,
00271 const EcalRecHitCollection* ecalEndcapHits,
00272 const edm::Handle<CaloTowerCollection> & hcalTowersHandle,
00273
00274 reco::VertexCollection & vertexCollection,
00275 reco::PhotonCollection & outputPhotonCollection, int& iSC,
00276 const EcalSeverityLevelAlgo * sevLv) {
00277
00278 const CaloGeometry* geometry = theCaloGeom_.product();
00279 const EcalRecHitCollection* hits = 0 ;
00280 std::vector<double> preselCutValues;
00281 std::vector<int> flags_, severitiesexcl_;
00282
00283 for(unsigned int lSC=0; lSC < photonCoreHandle->size(); lSC++) {
00284
00285 reco::PhotonCoreRef coreRef(reco::PhotonCoreRef(photonCoreHandle, lSC));
00286 reco::SuperClusterRef scRef=coreRef->superCluster();
00287
00288
00289
00290
00291 iSC++;
00292
00293 int subdet = scRef->seed()->hitsAndFractions()[0].first.subdetId();
00294 if (subdet==EcalBarrel) {
00295 preselCutValues = preselCutValuesBarrel_;
00296 hits = ecalBarrelHits;
00297 flags_ = flagsexclEB_;
00298 severitiesexcl_ = severitiesexclEB_;
00299 } else if (subdet==EcalEndcap) {
00300 preselCutValues = preselCutValuesEndcap_;
00301 hits = ecalEndcapHits;
00302 flags_ = flagsexclEE_;
00303 severitiesexcl_ = severitiesexclEE_;
00304 } else {
00305 edm::LogWarning("")<<"GEDPhotonProducer: do not know if it is a barrel or endcap SuperCluster";
00306 }
00307
00308
00309
00310 if (scRef->energy()/cosh(scRef->eta()) <= preselCutValues[0] ) continue;
00311
00312
00313 const CaloTowerCollection* hcalTowersColl = hcalTowersHandle.product();
00314 EgammaTowerIsolation towerIso1(hOverEConeSize_,0.,0.,1,hcalTowersColl) ;
00315 EgammaTowerIsolation towerIso2(hOverEConeSize_,0.,0.,2,hcalTowersColl) ;
00316 double HoE1=towerIso1.getTowerESum(&(*scRef))/scRef->energy();
00317 double HoE2=towerIso2.getTowerESum(&(*scRef))/scRef->energy();
00318
00319 EgammaHadTower towerIsoBehindClus(es);
00320 towerIsoBehindClus.setTowerCollection(hcalTowersHandle.product());
00321 std::vector<CaloTowerDetId> TowersBehindClus = towerIsoBehindClus.towersOf(*scRef);
00322 float hcalDepth1OverEcalBc = towerIsoBehindClus.getDepth1HcalESum(TowersBehindClus)/scRef->energy();
00323 float hcalDepth2OverEcalBc = towerIsoBehindClus.getDepth2HcalESum(TowersBehindClus)/scRef->energy();
00324
00325
00326
00327 float maxXtal = EcalClusterTools::eMax( *(scRef->seed()), &(*hits) );
00328
00329
00330 float e1x5 = EcalClusterTools::e1x5( *(scRef->seed()), &(*hits), &(*topology), flags_, severitiesexcl_, sevLv);
00331 float e2x5 = EcalClusterTools::e2x5Max( *(scRef->seed()), &(*hits), &(*topology),flags_, severitiesexcl_, sevLv );
00332 float e3x3 = EcalClusterTools::e3x3( *(scRef->seed()), &(*hits), &(*topology), flags_, severitiesexcl_, sevLv);
00333 float e5x5 = EcalClusterTools::e5x5( *(scRef->seed()), &(*hits), &(*topology),flags_, severitiesexcl_, sevLv);
00334 std::vector<float> cov = EcalClusterTools::covariances( *(scRef->seed()), &(*hits), &(*topology), geometry,flags_, severitiesexcl_, sevLv);
00335 std::vector<float> locCov = EcalClusterTools::localCovariances( *(scRef->seed()), &(*hits), &(*topology),flags_, severitiesexcl_, sevLv);
00336
00337 float sigmaEtaEta = sqrt(cov[0]);
00338 float sigmaIetaIeta = sqrt(locCov[0]);
00339
00340 math::XYZPoint caloPosition = scRef->position();
00341
00342
00344 double photonEnergy=1.;
00345 math::XYZPoint vtx(0.,0.,0.);
00346 if (vertexCollection.size()>0) vtx = vertexCollection.begin()->position();
00347
00348 math::XYZVector direction = caloPosition - vtx;
00349
00350 math::XYZVector momentum = direction.unit() ;
00351
00352
00353 math::XYZTLorentzVectorD p4(momentum.x(), momentum.y(), momentum.z(), photonEnergy );
00354 reco::Photon newCandidate(p4, caloPosition, coreRef, vtx);
00355
00356
00357
00358
00359
00360 reco::Photon::FiducialFlags fiducialFlags;
00361 reco::Photon::IsolationVariables isolVarR03, isolVarR04;
00362 thePhotonIsolationCalculator_-> calculate ( &newCandidate,evt,es,fiducialFlags,isolVarR04, isolVarR03);
00363 newCandidate.setFiducialVolumeFlags( fiducialFlags );
00364 newCandidate.setIsolationVariables(isolVarR04, isolVarR03 );
00365
00366
00367 reco::Photon::PflowIsolationVariables pfIso;
00368 reco::Photon::PflowIDVariables pfID;
00369 newCandidate.setPflowIsolationVariables(pfIso);
00370 newCandidate.setPflowIDVariables(pfID);
00371
00372
00374 reco::Photon::ShowerShape showerShape;
00375 showerShape.e1x5= e1x5;
00376 showerShape.e2x5= e2x5;
00377 showerShape.e3x3= e3x3;
00378 showerShape.e5x5= e5x5;
00379 showerShape.maxEnergyXtal = maxXtal;
00380 showerShape.sigmaEtaEta = sigmaEtaEta;
00381 showerShape.sigmaIetaIeta = sigmaIetaIeta;
00382 showerShape.hcalDepth1OverEcal = HoE1;
00383 showerShape.hcalDepth2OverEcal = HoE2;
00384 showerShape.hcalDepth1OverEcalBc = hcalDepth1OverEcalBc;
00385 showerShape.hcalDepth2OverEcalBc = hcalDepth2OverEcalBc;
00386 showerShape.hcalTowersBehindClusters = TowersBehindClus;
00387 newCandidate.setShowerShapeVariables ( showerShape );
00388
00391
00392 thePhotonEnergyCorrector_->calculate(evt, newCandidate, subdet, vertexCollection,es);
00393 if ( candidateP4type_ == "fromEcalEnergy") {
00394 newCandidate.setP4( newCandidate.p4(reco::Photon::ecal_photons) );
00395 newCandidate.setCandidateP4type(reco::Photon::ecal_photons);
00396 } else if ( candidateP4type_ == "fromRegression1") {
00397 newCandidate.setP4( newCandidate.p4(reco::Photon::regression1) );
00398 newCandidate.setCandidateP4type(reco::Photon::regression1);
00399 } else if ( candidateP4type_ == "fromRegression2") {
00400 newCandidate.setP4( newCandidate.p4(reco::Photon::regression2) );
00401 newCandidate.setCandidateP4type(reco::Photon::regression2);
00402 }
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412 reco::Photon::MIPVariables mipVar ;
00413 if(subdet==EcalBarrel && runMIPTagger_ )
00414 {
00415
00416 thePhotonMIPHaloTagger_-> MIPcalculate( &newCandidate,evt,es,mipVar);
00417 newCandidate.setMIPVariables(mipVar);
00418 }
00419
00420
00421
00423 bool isLooseEM=true;
00424 if ( newCandidate.pt() < highEt_) {
00425 if ( newCandidate.hadronicOverEm() >= preselCutValues[1] ) isLooseEM=false;
00426 if ( newCandidate.ecalRecHitSumEtConeDR04() > preselCutValues[2]+ preselCutValues[3]*newCandidate.pt() ) isLooseEM=false;
00427 if ( newCandidate.hcalTowerSumEtConeDR04() > preselCutValues[4]+ preselCutValues[5]*newCandidate.pt() ) isLooseEM=false;
00428 if ( newCandidate.nTrkSolidConeDR04() > int(preselCutValues[6]) ) isLooseEM=false;
00429 if ( newCandidate.nTrkHollowConeDR04() > int(preselCutValues[7]) ) isLooseEM=false;
00430 if ( newCandidate.trkSumPtSolidConeDR04() > preselCutValues[8] ) isLooseEM=false;
00431 if ( newCandidate.trkSumPtHollowConeDR04() > preselCutValues[9] ) isLooseEM=false;
00432 if ( newCandidate.sigmaIetaIeta() > preselCutValues[10] ) isLooseEM=false;
00433 }
00434
00435
00436
00437 if ( isLooseEM)
00438 outputPhotonCollection.push_back(newCandidate);
00439
00440
00441 }
00442 }
00443