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/PhotonProducer.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 PhotonProducer::PhotonProducer(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 PhotonProducer::~PhotonProducer()
00136 {
00137
00138
00139 }
00140
00141
00142
00143 void PhotonProducer::beginRun (edm::Run& 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 PhotonProducer::endRun (edm::Run& r, edm::EventSetup const & theEventSetup) {
00157
00158 delete thePhotonIsolationCalculator_;
00159 delete thePhotonMIPHaloTagger_;
00160 delete thePhotonEnergyCorrector_;
00161 }
00162
00163
00164 void PhotonProducer::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("PhotonProducer") << "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("PhotonProducer") << "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("PhotonProducer") << "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("PhotonProducer") << "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("PhotonProducer") << " 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 PhotonProducer::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 CaloSubdetectorGeometry* subDetGeometry =0 ;
00280 const CaloSubdetectorGeometry* geometryES = theCaloGeom_->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
00281 const EcalRecHitCollection* hits = 0 ;
00282 std::vector<double> preselCutValues;
00283 float minR9=0;
00284
00285
00286 std::vector<int> flags_, severitiesexcl_;
00287
00288 for(unsigned int lSC=0; lSC < photonCoreHandle->size(); lSC++) {
00289
00290 reco::PhotonCoreRef coreRef(reco::PhotonCoreRef(photonCoreHandle, lSC));
00291 reco::SuperClusterRef scRef=coreRef->superCluster();
00292
00293 iSC++;
00294
00295 int subdet = scRef->seed()->hitsAndFractions()[0].first.subdetId();
00296 subDetGeometry = theCaloGeom_->getSubdetectorGeometry(DetId::Ecal, subdet);
00297
00298 if (subdet==EcalBarrel) {
00299 preselCutValues = preselCutValuesBarrel_;
00300 minR9 = minR9Barrel_;
00301 hits = ecalBarrelHits;
00302 flags_ = flagsexclEB_;
00303 severitiesexcl_ = severitiesexclEB_;
00304 } else if (subdet==EcalEndcap) {
00305 preselCutValues = preselCutValuesEndcap_;
00306 minR9 = minR9Endcap_;
00307 hits = ecalEndcapHits;
00308 flags_ = flagsexclEE_;
00309 severitiesexcl_ = severitiesexclEE_;
00310 } else {
00311 edm::LogWarning("")<<"PhotonProducer: do not know if it is a barrel or endcap SuperCluster";
00312 }
00313
00314
00315
00316 if (scRef->energy()/cosh(scRef->eta()) <= preselCutValues[0] ) continue;
00317
00318
00319 const CaloTowerCollection* hcalTowersColl = hcalTowersHandle.product();
00320 EgammaTowerIsolation towerIso1(hOverEConeSize_,0.,0.,1,hcalTowersColl) ;
00321 EgammaTowerIsolation towerIso2(hOverEConeSize_,0.,0.,2,hcalTowersColl) ;
00322 double HoE1=towerIso1.getTowerESum(&(*scRef))/scRef->energy();
00323 double HoE2=towerIso2.getTowerESum(&(*scRef))/scRef->energy();
00324
00325 EgammaHadTower towerIsoBehindClus(es);
00326 towerIsoBehindClus.setTowerCollection(hcalTowersHandle.product());
00327 std::vector<CaloTowerDetId> TowersBehindClus = towerIsoBehindClus.towersOf(*scRef);
00328 float hcalDepth1OverEcalBc = towerIsoBehindClus.getDepth1HcalESum(TowersBehindClus)/scRef->energy();
00329 float hcalDepth2OverEcalBc = towerIsoBehindClus.getDepth2HcalESum(TowersBehindClus)/scRef->energy();
00330
00331
00332
00333
00334
00335 math::XYZPoint unconvPos = posCalculator_.Calculate_Location(scRef->seed()->hitsAndFractions(),hits,subDetGeometry,geometryES);
00336
00337
00338 float maxXtal = EcalClusterTools::eMax( *(scRef->seed()), &(*hits) );
00339
00340
00341 float e1x5 = EcalClusterTools::e1x5( *(scRef->seed()), &(*hits), &(*topology), flags_, severitiesexcl_, sevLv);
00342 float e2x5 = EcalClusterTools::e2x5Max( *(scRef->seed()), &(*hits), &(*topology),flags_, severitiesexcl_, sevLv );
00343 float e3x3 = EcalClusterTools::e3x3( *(scRef->seed()), &(*hits), &(*topology), flags_, severitiesexcl_, sevLv);
00344 float e5x5 = EcalClusterTools::e5x5( *(scRef->seed()), &(*hits), &(*topology),flags_, severitiesexcl_, sevLv);
00345 std::vector<float> cov = EcalClusterTools::covariances( *(scRef->seed()), &(*hits), &(*topology), geometry,flags_, severitiesexcl_, sevLv);
00346 std::vector<float> locCov = EcalClusterTools::localCovariances( *(scRef->seed()), &(*hits), &(*topology),flags_, severitiesexcl_, sevLv);
00347
00348 float sigmaEtaEta = sqrt(cov[0]);
00349 float sigmaIetaIeta = sqrt(locCov[0]);
00350 float r9 =e3x3/(scRef->rawEnergy());
00351
00352
00353 math::XYZPoint caloPosition;
00354 if (r9>minR9) {
00355 caloPosition = unconvPos;
00356 } else {
00357 caloPosition = scRef->position();
00358 }
00359
00361 double photonEnergy=1.;
00362 math::XYZPoint vtx(0.,0.,0.);
00363 if (vertexCollection.size()>0) vtx = vertexCollection.begin()->position();
00364
00365 math::XYZVector direction = caloPosition - vtx;
00366
00367 math::XYZVector momentum = direction.unit() ;
00368
00369
00370 math::XYZTLorentzVectorD p4(momentum.x(), momentum.y(), momentum.z(), photonEnergy );
00371 reco::Photon newCandidate(p4, caloPosition, coreRef, vtx);
00372
00373
00374
00375
00376 reco::Photon::FiducialFlags fiducialFlags;
00377 reco::Photon::IsolationVariables isolVarR03, isolVarR04;
00378 thePhotonIsolationCalculator_-> calculate ( &newCandidate,evt,es,fiducialFlags,isolVarR04, isolVarR03);
00379 newCandidate.setFiducialVolumeFlags( fiducialFlags );
00380 newCandidate.setIsolationVariables(isolVarR04, isolVarR03 );
00381
00382
00384 reco::Photon::ShowerShape showerShape;
00385 showerShape.e1x5= e1x5;
00386 showerShape.e2x5= e2x5;
00387 showerShape.e3x3= e3x3;
00388 showerShape.e5x5= e5x5;
00389 showerShape.maxEnergyXtal = maxXtal;
00390 showerShape.sigmaEtaEta = sigmaEtaEta;
00391 showerShape.sigmaIetaIeta = sigmaIetaIeta;
00392 showerShape.hcalDepth1OverEcal = HoE1;
00393 showerShape.hcalDepth2OverEcal = HoE2;
00394 showerShape.hcalDepth1OverEcalBc = hcalDepth1OverEcalBc;
00395 showerShape.hcalDepth2OverEcalBc = hcalDepth2OverEcalBc;
00396 showerShape.hcalTowersBehindClusters = TowersBehindClus;
00397 newCandidate.setShowerShapeVariables ( showerShape );
00398
00401
00402 thePhotonEnergyCorrector_->calculate(evt, newCandidate, subdet, vertexCollection,es);
00403 if ( candidateP4type_ == "fromEcalEnergy") {
00404 newCandidate.setP4( newCandidate.p4(reco::Photon::ecal_photons) );
00405 newCandidate.setCandidateP4type(reco::Photon::ecal_photons);
00406 } else if ( candidateP4type_ == "fromRegression") {
00407 newCandidate.setP4( newCandidate.p4(reco::Photon::regression1) );
00408 newCandidate.setCandidateP4type(reco::Photon::regression1);
00409 }
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419 reco::Photon::MIPVariables mipVar ;
00420 if(subdet==EcalBarrel && runMIPTagger_ )
00421 {
00422
00423 thePhotonMIPHaloTagger_-> MIPcalculate( &newCandidate,evt,es,mipVar);
00424 newCandidate.setMIPVariables(mipVar);
00425 }
00426
00427
00428
00430 bool isLooseEM=true;
00431 if ( newCandidate.pt() < highEt_) {
00432 if ( newCandidate.hadronicOverEm() >= preselCutValues[1] ) isLooseEM=false;
00433 if ( newCandidate.ecalRecHitSumEtConeDR04() > preselCutValues[2]+ preselCutValues[3]*newCandidate.pt() ) isLooseEM=false;
00434 if ( newCandidate.hcalTowerSumEtConeDR04() > preselCutValues[4]+ preselCutValues[5]*newCandidate.pt() ) isLooseEM=false;
00435 if ( newCandidate.nTrkSolidConeDR04() > int(preselCutValues[6]) ) isLooseEM=false;
00436 if ( newCandidate.nTrkHollowConeDR04() > int(preselCutValues[7]) ) isLooseEM=false;
00437 if ( newCandidate.trkSumPtSolidConeDR04() > preselCutValues[8] ) isLooseEM=false;
00438 if ( newCandidate.trkSumPtHollowConeDR04() > preselCutValues[9] ) isLooseEM=false;
00439 if ( newCandidate.sigmaIetaIeta() > preselCutValues[10] ) isLooseEM=false;
00440 }
00441
00442
00443
00444 if ( isLooseEM)
00445 outputPhotonCollection.push_back(newCandidate);
00446
00447
00448 }
00449 }
00450