Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-03-03 23:53:12 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!