EgammaPhotonProducers/src/PhotonProducer.cc

0001 /** \class PhotonProducer
0002  **
0003  **
0004  **  \author Nancy Marinelli, U. of Notre Dame, US
0005  **
0006  ***/
0007
0008 #include "CommonTools/Utils/interface/StringToEnumValue.h"
0009 #include "CondFormats/EcalObjects/interface/EcalFunctionParameters.h"
0010 #include "DataFormats/Common/interface/Handle.h"
0011 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
0012 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0013 #include "DataFormats/EgammaCandidates/interface/PhotonCore.h"
0014 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0015 #include "DataFormats/EgammaReco/interface/BasicClusterShapeAssociation.h"
0016 #include "DataFormats/EgammaReco/interface/ClusterShape.h"
0017 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0018 #include "DataFormats/VertexReco/interface/Vertex.h"
0019 #include "FWCore/Framework/interface/ESHandle.h"
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/Framework/interface/stream/EDProducer.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/Utilities/interface/Exception.h"
0026 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0027 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0028 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0029 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0030 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0031 #include "RecoCaloTools/Selectors/interface/CaloConeSelector.h"
0032 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0033 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
0034 #include "RecoEgamma/EgammaPhotonAlgos/interface/PhotonEnergyCorrector.h"
0035 #include "RecoEgamma/PhotonIdentification/interface/PhotonIsolationCalculator.h"
0036 #include "RecoEgamma/PhotonIdentification/interface/PhotonMIPHaloTagger.h"
0037 #include "RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h"
0038 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0039 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0040 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronHcalHelper.h"
0041 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0042 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0043 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0044
0045 #include <vector>
0046
0047 // PhotonProducer inherits from EDProducer, so it can be a module:
0048 class PhotonProducer : public edm::stream::EDProducer<> {
0049 public:
0050   PhotonProducer(const edm::ParameterSet& ps);
0051
0052   void produce(edm::Event& evt, const edm::EventSetup& es) override;
0053
0054 private:
0055   void fillPhotonCollection(edm::Event& evt,
0056                             edm::EventSetup const& es,
0057                             const edm::Handle<reco::PhotonCoreCollection>& photonCoreHandle,
0058                             const CaloTopology* topology,
0059                             const HcalPFCuts* hcalCuts,
0060                             const EcalRecHitCollection* ecalBarrelHits,
0061                             const EcalRecHitCollection* ecalEndcapHits,
0062                             ElectronHcalHelper const& hcalHelperCone,
0063                             ElectronHcalHelper const& hcalHelperBc,
0064                             reco::VertexCollection& pvVertices,
0065                             reco::PhotonCollection& outputCollection,
0066                             int& iSC);
0067
0068   // std::string PhotonCoreCollection_;
0069   std::string PhotonCollection_;
0070   edm::EDGetTokenT<reco::PhotonCoreCollection> photonCoreProducer_;
0071   edm::EDGetTokenT<EcalRecHitCollection> barrelEcalHits_;
0072   edm::EDGetTokenT<EcalRecHitCollection> endcapEcalHits_;
0073   edm::EDGetTokenT<HBHERecHitCollection> hbheRecHits_;
0074   edm::EDGetTokenT<reco::VertexCollection> vertexProducer_;
0075   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0076   const edm::ESGetToken<CaloTopology, CaloTopologyRecord> topologyToken_;
0077
0078   //AA
0079   //Flags and severities to be excluded from calculations
0080
0081   std::vector<int> flagsexclEB_;
0082   std::vector<int> flagsexclEE_;
0083   std::vector<int> severitiesexclEB_;
0084   std::vector<int> severitiesexclEE_;
0085
0086   double hOverEConeSize_;
0087   double highEt_;
0088   double minR9Barrel_;
0089   double minR9Endcap_;
0090   bool runMIPTagger_;
0091   bool runMVABasedHaloTagger_;
0092
0093   bool validConversions_;
0094
0095   bool usePrimaryVertex_;
0096
0097   PositionCalc posCalculator_;
0098
0099   bool validPixelSeeds_;
0100   PhotonIsolationCalculator photonIsolationCalculator_;
0101
0102   //MIP
0103   const PhotonMIPHaloTagger photonMIPHaloTagger_;
0104   //MVA based Halo tagger for the EE photons
0105   std::unique_ptr<const PhotonMVABasedHaloTagger> photonMVABasedHaloTagger_ = nullptr;
0106
0107   std::vector<double> preselCutValuesBarrel_;
0108   std::vector<double> preselCutValuesEndcap_;
0109
0110   PhotonEnergyCorrector photonEnergyCorrector_;
0111   std::string candidateP4type_;
0112
0113   // additional configuration and helpers
0114   std::unique_ptr<ElectronHcalHelper> hcalHelperCone_;
0115   std::unique_ptr<ElectronHcalHelper> hcalHelperBc_;
0116   bool hcalRun2EffDepth_;
0117
0118   edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd> hcalCutsToken_;
0119   bool cutsFromDB_;
0120 };
0121
0122 #include "FWCore/Framework/interface/MakerMacros.h"
0123 DEFINE_FWK_MODULE(PhotonProducer);
0124
0125 PhotonProducer::PhotonProducer(const edm::ParameterSet& config)
0126     : caloGeomToken_(esConsumes()),
0127       topologyToken_(esConsumes()),
0128       flagsexclEB_{StringToEnumValue<EcalRecHit::Flags>(
0129           config.getParameter<std::vector<std::string>>("RecHitFlagToBeExcludedEB"))},
0130       flagsexclEE_{StringToEnumValue<EcalRecHit::Flags>(
0131           config.getParameter<std::vector<std::string>>("RecHitFlagToBeExcludedEE"))},
0132       severitiesexclEB_{StringToEnumValue<EcalSeverityLevel::SeverityLevel>(
0133           config.getParameter<std::vector<std::string>>("RecHitSeverityToBeExcludedEB"))},
0134       severitiesexclEE_{StringToEnumValue<EcalSeverityLevel::SeverityLevel>(
0135           config.getParameter<std::vector<std::string>>("RecHitSeverityToBeExcludedEE"))},
0136       photonIsolationCalculator_(config.getParameter<edm::ParameterSet>("isolationSumsCalculatorSet"),
0137                                  flagsexclEB_,
0138                                  flagsexclEE_,
0139                                  severitiesexclEB_,
0140                                  severitiesexclEE_,
0141                                  consumesCollector()),
0142       photonMIPHaloTagger_(config.getParameter<edm::ParameterSet>("mipVariableSet"), consumesCollector()),
0143       photonEnergyCorrector_(config, consumesCollector()) {
0144   // use configuration file to setup input/output collection names
0145
0146   photonCoreProducer_ = consumes<reco::PhotonCoreCollection>(config.getParameter<edm::InputTag>("photonCoreProducer"));
0147   barrelEcalHits_ = consumes<EcalRecHitCollection>(config.getParameter<edm::InputTag>("barrelEcalHits"));
0148   endcapEcalHits_ = consumes<EcalRecHitCollection>(config.getParameter<edm::InputTag>("endcapEcalHits"));
0149   vertexProducer_ = consumes<reco::VertexCollection>(config.getParameter<edm::InputTag>("primaryVertexProducer"));
0150   hbheRecHits_ = consumes<HBHERecHitCollection>(config.getParameter<edm::InputTag>("hbheRecHits"));
0151   hOverEConeSize_ = config.getParameter<double>("hOverEConeSize");
0152   highEt_ = config.getParameter<double>("highEt");
0153   // R9 value to decide converted/unconverted
0154   minR9Barrel_ = config.getParameter<double>("minR9Barrel");
0155   minR9Endcap_ = config.getParameter<double>("minR9Endcap");
0156   usePrimaryVertex_ = config.getParameter<bool>("usePrimaryVertex");
0157   runMIPTagger_ = config.getParameter<bool>("runMIPTagger");
0158   runMVABasedHaloTagger_ = config.getParameter<bool>("runMVABasedHaloTagger");
0159
0160   candidateP4type_ = config.getParameter<std::string>("candidateP4type");
0161
0162   edm::ParameterSet posCalcParameters = config.getParameter<edm::ParameterSet>("posCalcParameters");
0163   posCalculator_ = PositionCalc(posCalcParameters);
0164
0165   //Retrieve HCAL PF thresholds - from config or from DB
0166   cutsFromDB_ = config.getParameter<bool>("usePFThresholdsFromDB");
0167   if (cutsFromDB_) {
0168     hcalCutsToken_ = esConsumes<HcalPFCuts, HcalPFCutsRcd>(edm::ESInputTag("", "withTopo"));
0169   }
0170
0171   ElectronHcalHelper::Configuration cfgCone, cfgBc;
0172   cfgCone.hOverEConeSize = hOverEConeSize_;
0173   if (cfgCone.hOverEConeSize > 0) {
0174     cfgCone.onlyBehindCluster = false;
0175     cfgCone.checkHcalStatus = false;
0176
0177     cfgCone.hbheRecHits = hbheRecHits_;
0178
0179     cfgCone.eThresHB = config.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
0180     cfgCone.maxSeverityHB = config.getParameter<int>("maxHcalRecHitSeverity");
0181     cfgCone.eThresHE = config.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
0182     cfgCone.maxSeverityHE = cfgCone.maxSeverityHB;
0183   }
0184
0185   cfgBc.hOverEConeSize = 0.;
0186   cfgBc.onlyBehindCluster = true;
0187   cfgBc.checkHcalStatus = false;
0188
0189   cfgBc.hbheRecHits = hbheRecHits_;
0190
0191   cfgBc.eThresHB = config.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
0192   cfgBc.maxSeverityHB = config.getParameter<int>("maxHcalRecHitSeverity");
0193   cfgBc.eThresHE = config.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
0194   cfgBc.maxSeverityHE = cfgBc.maxSeverityHB;
0195
0196   hcalHelperCone_ = std::make_unique<ElectronHcalHelper>(cfgCone, consumesCollector());
0197   hcalHelperBc_ = std::make_unique<ElectronHcalHelper>(cfgBc, consumesCollector());
0198
0199   hcalRun2EffDepth_ = config.getParameter<bool>("hcalRun2EffDepth");
0200
0201   //AA
0202
0203   //
0204
0205   // Parameters for the position calculation:
0206   //  std::map<std::string,double> providedParameters;
0207   // providedParameters.insert(std::make_pair("LogWeighted",config.getParameter<bool>("posCalc_logweight")));
0208   //providedParameters.insert(std::make_pair("T0_barl",config.getParameter<double>("posCalc_t0_barl")));
0209   //providedParameters.insert(std::make_pair("T0_endc",config.getParameter<double>("posCalc_t0_endc")));
0210   //providedParameters.insert(std::make_pair("T0_endcPresh",config.getParameter<double>("posCalc_t0_endcPresh")));
0211   //providedParameters.insert(std::make_pair("W0",config.getParameter<double>("posCalc_w0")));
0212   //providedParameters.insert(std::make_pair("X0",config.getParameter<double>("posCalc_x0")));
0213   //posCalculator_ = PositionCalc(providedParameters);
0214   // cut values for pre-selection
0215   preselCutValuesBarrel_.push_back(config.getParameter<double>("minSCEtBarrel"));
0216   preselCutValuesBarrel_.push_back(config.getParameter<double>("maxHoverEBarrel"));
0217   preselCutValuesBarrel_.push_back(config.getParameter<double>("ecalRecHitSumEtOffsetBarrel"));
0218   preselCutValuesBarrel_.push_back(config.getParameter<double>("ecalRecHitSumEtSlopeBarrel"));
0219   preselCutValuesBarrel_.push_back(config.getParameter<double>("hcalRecHitSumEtOffsetBarrel"));
0220   preselCutValuesBarrel_.push_back(config.getParameter<double>("hcalRecHitSumEtSlopeBarrel"));
0221   preselCutValuesBarrel_.push_back(config.getParameter<double>("nTrackSolidConeBarrel"));
0222   preselCutValuesBarrel_.push_back(config.getParameter<double>("nTrackHollowConeBarrel"));
0223   preselCutValuesBarrel_.push_back(config.getParameter<double>("trackPtSumSolidConeBarrel"));
0224   preselCutValuesBarrel_.push_back(config.getParameter<double>("trackPtSumHollowConeBarrel"));
0225   preselCutValuesBarrel_.push_back(config.getParameter<double>("sigmaIetaIetaCutBarrel"));
0226   //
0227   preselCutValuesEndcap_.push_back(config.getParameter<double>("minSCEtEndcap"));
0228   preselCutValuesEndcap_.push_back(config.getParameter<double>("maxHoverEEndcap"));
0229   preselCutValuesEndcap_.push_back(config.getParameter<double>("ecalRecHitSumEtOffsetEndcap"));
0230   preselCutValuesEndcap_.push_back(config.getParameter<double>("ecalRecHitSumEtSlopeEndcap"));
0231   preselCutValuesEndcap_.push_back(config.getParameter<double>("hcalRecHitSumEtOffsetEndcap"));
0232   preselCutValuesEndcap_.push_back(config.getParameter<double>("hcalRecHitSumEtSlopeEndcap"));
0233   preselCutValuesEndcap_.push_back(config.getParameter<double>("nTrackSolidConeEndcap"));
0234   preselCutValuesEndcap_.push_back(config.getParameter<double>("nTrackHollowConeEndcap"));
0235   preselCutValuesEndcap_.push_back(config.getParameter<double>("trackPtSumSolidConeEndcap"));
0236   preselCutValuesEndcap_.push_back(config.getParameter<double>("trackPtSumHollowConeEndcap"));
0237   preselCutValuesEndcap_.push_back(config.getParameter<double>("sigmaIetaIetaCutEndcap"));
0238   //
0239
0240   // Register the product
0241   produces<reco::PhotonCollection>(PhotonCollection_);
0242 }
0243
0244 void PhotonProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
0245   HcalPFCuts const* hcalCuts = nullptr;
0246   if (cutsFromDB_) {
0247     hcalCuts = &theEventSetup.getData(hcalCutsToken_);
0248   }
0249   using namespace edm;
0250   //  nEvt_++;
0251
0252   reco::PhotonCollection outputPhotonCollection;
0253   auto outputPhotonCollection_p = std::make_unique<reco::PhotonCollection>();
0254
0255   // Get the PhotonCore collection
0256   bool validPhotonCoreHandle = true;
0257   Handle<reco::PhotonCoreCollection> photonCoreHandle;
0258   theEvent.getByToken(photonCoreProducer_, photonCoreHandle);
0259   if (!photonCoreHandle.isValid()) {
0260     edm::LogError("PhotonProducer") << "Error! Can't get the photonCoreProducer";
0261     validPhotonCoreHandle = false;
0262   }
0263
0264   // Get EcalRecHits
0265   bool validEcalRecHits = true;
0266   Handle<EcalRecHitCollection> barrelHitHandle;
0267   EcalRecHitCollection barrelRecHits;
0268   theEvent.getByToken(barrelEcalHits_, barrelHitHandle);
0269   if (!barrelHitHandle.isValid()) {
0270     edm::LogError("PhotonProducer") << "Error! Can't get the barrelEcalHits";
0271     validEcalRecHits = false;
0272   }
0273   if (validEcalRecHits)
0274     barrelRecHits = *(barrelHitHandle.product());
0275
0276   Handle<EcalRecHitCollection> endcapHitHandle;
0277   theEvent.getByToken(endcapEcalHits_, endcapHitHandle);
0278   EcalRecHitCollection endcapRecHits;
0279   if (!endcapHitHandle.isValid()) {
0280     edm::LogError("PhotonProducer") << "Error! Can't get the endcapEcalHits";
0281     validEcalRecHits = false;
0282   }
0283   if (validEcalRecHits)
0284     endcapRecHits = *(endcapHitHandle.product());
0285
0286   const CaloTopology* topology = &theEventSetup.getData(topologyToken_);
0287
0288   // prepare access to hcal data
0289   hcalHelperCone_->beginEvent(theEvent, theEventSetup);
0290   hcalHelperBc_->beginEvent(theEvent, theEventSetup);
0291
0292   // Get the primary event vertex
0293   Handle<reco::VertexCollection> vertexHandle;
0294   reco::VertexCollection vertexCollection;
0295   bool validVertex = true;
0296   if (usePrimaryVertex_) {
0297     theEvent.getByToken(vertexProducer_, vertexHandle);
0298     if (!vertexHandle.isValid()) {
0299       edm::LogWarning("PhotonProducer") << "Error! Can't get the product primary Vertex Collection "
0300                                         << "\n";
0301       validVertex = false;
0302     }
0303     if (validVertex)
0304       vertexCollection = *(vertexHandle.product());
0305   }
0306
0307   int iSC = 0;  // index in photon collection
0308   // Loop over barrel and endcap SC collections and fill the  photon collection
0309   if (validPhotonCoreHandle)
0310     fillPhotonCollection(theEvent,
0311                          theEventSetup,
0312                          photonCoreHandle,
0313                          topology,
0314                          hcalCuts,
0315                          &barrelRecHits,
0316                          &endcapRecHits,
0317                          *hcalHelperCone_,
0318                          *hcalHelperBc_,
0319                          vertexCollection,
0320                          outputPhotonCollection,
0321                          iSC);
0322
0323   // put the product in the event
0324   edm::LogInfo("PhotonProducer") << " Put in the event " << iSC << " Photon Candidates \n";
0325   outputPhotonCollection_p->assign(outputPhotonCollection.begin(), outputPhotonCollection.end());
0326
0327   // go back to run2-like 2 effective depths if desired - depth 1 is the normal depth 1, depth 2 is the sum over the rest
0328   if (hcalRun2EffDepth_) {
0329     for (auto& pho : *outputPhotonCollection_p)
0330       pho.hcalToRun2EffDepth();
0331   }
0332
0333   theEvent.put(std::move(outputPhotonCollection_p), PhotonCollection_);
0334 }
0335
0336 void PhotonProducer::fillPhotonCollection(edm::Event& evt,
0337                                           edm::EventSetup const& es,
0338                                           const edm::Handle<reco::PhotonCoreCollection>& photonCoreHandle,
0339                                           const CaloTopology* topology,
0340                                           const HcalPFCuts* hcalCuts,
0341                                           const EcalRecHitCollection* ecalBarrelHits,
0342                                           const EcalRecHitCollection* ecalEndcapHits,
0343                                           ElectronHcalHelper const& hcalHelperCone,
0344                                           ElectronHcalHelper const& hcalHelperBc,
0345                                           reco::VertexCollection& vertexCollection,
0346                                           reco::PhotonCollection& outputPhotonCollection,
0347                                           int& iSC) {
0348   // get the geometry from the event setup:
0349   const CaloGeometry* geometry = &es.getData(caloGeomToken_);
0350   const CaloSubdetectorGeometry* subDetGeometry = nullptr;
0351   const CaloSubdetectorGeometry* geometryES = geometry->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0352   const EcalRecHitCollection* hits = nullptr;
0353   std::vector<double> preselCutValues;
0354   float minR9 = 0;
0355
0356   photonEnergyCorrector_.init(es);
0357
0358   std::vector<int> flags_, severitiesexcl_;
0359
0360   for (unsigned int lSC = 0; lSC < photonCoreHandle->size(); lSC++) {
0361     reco::PhotonCoreRef coreRef(reco::PhotonCoreRef(photonCoreHandle, lSC));
0362     reco::SuperClusterRef scRef = coreRef->superCluster();
0363     iSC++;
0364
0365     int subdet = scRef->seed()->hitsAndFractions()[0].first.subdetId();
0366     subDetGeometry = geometry->getSubdetectorGeometry(DetId::Ecal, subdet);
0367
0368     if (subdet == EcalBarrel) {
0369       preselCutValues = preselCutValuesBarrel_;
0370       minR9 = minR9Barrel_;
0371       hits = ecalBarrelHits;
0372       flags_ = flagsexclEB_;
0373       severitiesexcl_ = severitiesexclEB_;
0374     } else if (subdet == EcalEndcap) {
0375       preselCutValues = preselCutValuesEndcap_;
0376       minR9 = minR9Endcap_;
0377       hits = ecalEndcapHits;
0378       flags_ = flagsexclEE_;
0379       severitiesexcl_ = severitiesexclEE_;
0380     } else {
0381       edm::LogWarning("") << "PhotonProducer: do not know if it is a barrel or endcap SuperCluster";
0382     }
0383     if (hits == nullptr)
0384       continue;
0385
0386     // SC energy preselection
0387     if (scRef->energy() / cosh(scRef->eta()) <= preselCutValues[0])
0388       continue;
0389
0390     // recalculate position of seed BasicCluster taking shower depth for unconverted photon
0391     math::XYZPoint unconvPos =
0392         posCalculator_.Calculate_Location(scRef->seed()->hitsAndFractions(), hits, subDetGeometry, geometryES);
0393
0394     float maxXtal = EcalClusterTools::eMax(*(scRef->seed()), &(*hits));
0395     //AA
0396     //Change these to consider severity level of hits
0397     float e1x5 = EcalClusterTools::e1x5(*(scRef->seed()), &(*hits), &(*topology));
0398     float e2x5 = EcalClusterTools::e2x5Max(*(scRef->seed()), &(*hits), &(*topology));
0399     float e3x3 = EcalClusterTools::e3x3(*(scRef->seed()), &(*hits), &(*topology));
0400     float e5x5 = EcalClusterTools::e5x5(*(scRef->seed()), &(*hits), &(*topology));
0401     const auto& cov = EcalClusterTools::covariances(*(scRef->seed()), &(*hits), &(*topology), geometry);
0402     const auto& locCov = EcalClusterTools::localCovariances(*(scRef->seed()), &(*hits), &(*topology));
0403
0404     float sigmaEtaEta = sqrt(cov[0]);
0405     float sigmaIetaIeta = sqrt(locCov[0]);
0406     float r9 = e3x3 / (scRef->rawEnergy());
0407
0408     float full5x5_maxXtal = noZS::EcalClusterTools::eMax(*(scRef->seed()), &(*hits));
0409     //AA
0410     //Change these to consider severity level of hits
0411     float full5x5_e1x5 = noZS::EcalClusterTools::e1x5(*(scRef->seed()), &(*hits), &(*topology));
0412     float full5x5_e2x5 = noZS::EcalClusterTools::e2x5Max(*(scRef->seed()), &(*hits), &(*topology));
0413     float full5x5_e3x3 = noZS::EcalClusterTools::e3x3(*(scRef->seed()), &(*hits), &(*topology));
0414     float full5x5_e5x5 = noZS::EcalClusterTools::e5x5(*(scRef->seed()), &(*hits), &(*topology));
0415     const auto& full5x5_cov = noZS::EcalClusterTools::covariances(*(scRef->seed()), &(*hits), &(*topology), geometry);
0416     const auto& full5x5_locCov = noZS::EcalClusterTools::localCovariances(*(scRef->seed()), &(*hits), &(*topology));
0417
0418     float full5x5_sigmaEtaEta = sqrt(full5x5_cov[0]);
0419     float full5x5_sigmaIetaIeta = sqrt(full5x5_locCov[0]);
0420
0421     // compute position of ECAL shower
0422     math::XYZPoint caloPosition;
0423     if (r9 > minR9) {
0424       caloPosition = unconvPos;
0425     } else {
0426       caloPosition = scRef->position();
0427     }
0428
0429     //// energy determination -- Default to create the candidate. Afterwards corrections are applied
0430     double photonEnergy = 1.;
0431     math::XYZPoint vtx(0., 0., 0.);
0432     if (!vertexCollection.empty())
0433       vtx = vertexCollection.begin()->position();
0434     // compute momentum vector of photon from primary vertex and cluster position
0435     math::XYZVector direction = caloPosition - vtx;
0436     math::XYZVector momentum = direction.unit();
0437
0438     // Create dummy candidate with unit momentum and zero energy to allow setting of all variables. The energy is set for last.
0439     math::XYZTLorentzVectorD p4(momentum.x(), momentum.y(), momentum.z(), photonEnergy);
0440     reco::Photon newCandidate(p4, caloPosition, coreRef, vtx);
0441
0442     // Calculate fiducial flags and isolation variable. Blocked are filled from the isolationCalculator
0443     reco::Photon::FiducialFlags fiducialFlags;
0444     reco::Photon::IsolationVariables isolVarR03, isolVarR04;
0445     photonIsolationCalculator_.calculate(&newCandidate, evt, es, fiducialFlags, isolVarR04, isolVarR03, hcalCuts);
0446     newCandidate.setFiducialVolumeFlags(fiducialFlags);
0447     newCandidate.setIsolationVariables(isolVarR04, isolVarR03);
0448
0449     /// fill shower shape block
0450     reco::Photon::ShowerShape showerShape;
0451     showerShape.e1x5 = e1x5;
0452     showerShape.e2x5 = e2x5;
0453     showerShape.e3x3 = e3x3;
0454     showerShape.e5x5 = e5x5;
0455     showerShape.maxEnergyXtal = maxXtal;
0456     showerShape.sigmaEtaEta = sigmaEtaEta;
0457     showerShape.sigmaIetaIeta = sigmaIetaIeta;
0458     for (uint id = 0; id < showerShape.hcalOverEcal.size(); ++id) {
0459       showerShape.hcalOverEcal[id] = hcalHelperCone.hcalESum(*scRef, id + 1, hcalCuts) / scRef->energy();
0460       showerShape.hcalOverEcalBc[id] = hcalHelperBc.hcalESum(*scRef, id + 1, hcalCuts) / scRef->energy();
0461     }
0462     showerShape.hcalTowersBehindClusters = hcalHelperBc.hcalTowersBehindClusters(*scRef);
0463     showerShape.pre7DepthHcal = false;
0464     newCandidate.setShowerShapeVariables(showerShape);
0465
0466     /// fill full5x5 shower shape block
0467     reco::Photon::ShowerShape full5x5_showerShape;
0468     full5x5_showerShape.e1x5 = full5x5_e1x5;
0469     full5x5_showerShape.e2x5 = full5x5_e2x5;
0470     full5x5_showerShape.e3x3 = full5x5_e3x3;
0471     full5x5_showerShape.e5x5 = full5x5_e5x5;
0472     full5x5_showerShape.maxEnergyXtal = full5x5_maxXtal;
0473     full5x5_showerShape.sigmaEtaEta = full5x5_sigmaEtaEta;
0474     full5x5_showerShape.sigmaIetaIeta = full5x5_sigmaIetaIeta;
0475     for (uint id = 0; id < full5x5_showerShape.hcalOverEcal.size(); ++id) {
0476       full5x5_showerShape.hcalOverEcal[id] = hcalHelperCone.hcalESum(*scRef, id + 1, hcalCuts) / full5x5_e5x5;
0477       full5x5_showerShape.hcalOverEcalBc[id] = hcalHelperBc.hcalESum(*scRef, id + 1, hcalCuts) / full5x5_e5x5;
0478     }
0479     full5x5_showerShape.hcalTowersBehindClusters = hcalHelperBc.hcalTowersBehindClusters(*scRef);
0480     full5x5_showerShape.pre7DepthHcal = false;
0481     newCandidate.full5x5_setShowerShapeVariables(full5x5_showerShape);
0482
0483     /// get ecal photon specific corrected energy
0484     /// plus values from regressions     and store them in the Photon
0485     // Photon candidate takes by default (set in photons_cfi.py)  a 4-momentum derived from the ecal photon-specific corrections.
0486     photonEnergyCorrector_.calculate(evt, newCandidate, subdet, vertexCollection, es);
0487     if (candidateP4type_ == "fromEcalEnergy") {
0488       newCandidate.setP4(newCandidate.p4(reco::Photon::ecal_photons));
0489       newCandidate.setCandidateP4type(reco::Photon::ecal_photons);
0490     } else if (candidateP4type_ == "fromRegression") {
0491       newCandidate.setP4(newCandidate.p4(reco::Photon::regression1));
0492       newCandidate.setCandidateP4type(reco::Photon::regression1);
0493     }
0494
0495     // fill MIP Vairables for Halo: Block for MIP are filled from PhotonMIPHaloTagger
0496     if (subdet == EcalBarrel && runMIPTagger_) {
0497       auto mipVar = photonMIPHaloTagger_.mipCalculate(newCandidate, evt, es);
0498       newCandidate.setMIPVariables(mipVar);
0499     }
0500
0501     /// Pre-selection loose  isolation cuts
0502     bool isLooseEM = true;
0503     if (newCandidate.pt() < highEt_) {
0504       if (newCandidate.hadronicOverEm() >= preselCutValues[1])
0505         isLooseEM = false;
0506       if (newCandidate.ecalRecHitSumEtConeDR04() > preselCutValues[2] + preselCutValues[3] * newCandidate.pt())
0507         isLooseEM = false;
0508       if (newCandidate.hcalTowerSumEtConeDR04() > preselCutValues[4] + preselCutValues[5] * newCandidate.pt())
0509         isLooseEM = false;
0510       if (newCandidate.nTrkSolidConeDR04() > int(preselCutValues[6]))
0511         isLooseEM = false;
0512       if (newCandidate.nTrkHollowConeDR04() > int(preselCutValues[7]))
0513         isLooseEM = false;
0514       if (newCandidate.trkSumPtSolidConeDR04() > preselCutValues[8])
0515         isLooseEM = false;
0516       if (newCandidate.trkSumPtHollowConeDR04() > preselCutValues[9])
0517         isLooseEM = false;
0518       if (newCandidate.sigmaIetaIeta() > preselCutValues[10])
0519         isLooseEM = false;
0520     }
0521
0522     if (isLooseEM)
0523       outputPhotonCollection.push_back(newCandidate);
0524   }
0525 }