Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-07-15 08:47:39 

0001 #include "DQMOffline/Trigger/interface/EgHLTOffHelper.h"
0002 
0003 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0004 #include "DataFormats/DetId/interface/DetId.h"
0005 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0006 #include "FWCore/Common/interface/TriggerNames.h"
0007 #include "FWCore/Framework/interface/ConsumesCollector.h"
0008 
0009 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0010 
0011 #include "DataFormats/TrackReco/interface/Track.h"
0012 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0013 
0014 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHLTTrackIsolation.h"
0015 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTowerIsolation.h"
0016 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaRecHitIsolation.h"
0017 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0018 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0019 
0020 #include "DQMOffline/Trigger/interface/EgHLTTrigCodes.h"
0021 #include "DQMOffline/Trigger/interface/EgHLTTrigTools.h"
0022 #include "DQMOffline/Trigger/interface/EgHLTErrCodes.h"
0023 
0024 #include <iostream>
0025 
0026 using namespace egHLT;
0027 
0028 OffHelper::~OffHelper() {
0029   if (hltEleTrkIsolAlgo_)
0030     delete hltEleTrkIsolAlgo_;
0031   if (hltPhoTrkIsolAlgo_)
0032     delete hltPhoTrkIsolAlgo_;
0033 }
0034 
0035 void OffHelper::setup(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC) {
0036   ecalRecHitsEBToken = iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("BarrelRecHitCollection"));
0037   ecalRecHitsEEToken = iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("EndcapRecHitCollection"));
0038   caloJetsToken = iC.consumes<reco::CaloJetCollection>(conf.getParameter<edm::InputTag>("CaloJetCollection"));
0039   isolTrkToken = iC.consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("IsolTrackCollection"));
0040   hbheHitsToken = iC.consumes<HBHERecHitCollection>(conf.getParameter<edm::InputTag>("HBHERecHitCollection"));
0041   hfHitsToken = iC.consumes<HFRecHitCollection>(conf.getParameter<edm::InputTag>("HFRecHitCollection"));
0042   electronsToken = iC.consumes<reco::GsfElectronCollection>(conf.getParameter<edm::InputTag>("ElectronCollection"));
0043   photonsToken = iC.consumes<reco::PhotonCollection>(conf.getParameter<edm::InputTag>("PhotonCollection"));
0044   triggerSummaryToken = iC.consumes<trigger::TriggerEvent>(conf.getParameter<edm::InputTag>("triggerSummaryLabel"));
0045   hltTag_ = conf.getParameter<std::string>("hltTag");
0046   beamSpotToken = iC.consumes<reco::BeamSpot>(conf.getParameter<edm::InputTag>("BeamSpotProducer"));
0047   caloTowersToken = iC.consumes<CaloTowerCollection>(conf.getParameter<edm::InputTag>("CaloTowers"));
0048   trigResultsToken = iC.consumes<edm::TriggerResults>(conf.getParameter<edm::InputTag>("TrigResults"));
0049   vertexToken = iC.consumes<reco::VertexCollection>(conf.getParameter<edm::InputTag>("VertexCollection"));
0050 
0051   caloGeomToken_ = iC.esConsumes();
0052   caloTopoToken_ = iC.esConsumes();
0053   magFieldToken_ = iC.esConsumes();
0054   ecalSeverityToken_ = iC.esConsumes();
0055 
0056   eleCuts_.setup(conf.getParameter<edm::ParameterSet>("eleCuts"));
0057   eleLooseCuts_.setup(conf.getParameter<edm::ParameterSet>("eleLooseCuts"));
0058   phoCuts_.setup(conf.getParameter<edm::ParameterSet>("phoCuts"));
0059   phoLooseCuts_.setup(conf.getParameter<edm::ParameterSet>("phoLooseCuts"));
0060 
0061   //now we have the isolations completely configurable via python
0062   hltEMIsolOuterCone_ = conf.getParameter<double>("hltEMIsolOuterCone");
0063   hltEMIsolInnerConeEB_ = conf.getParameter<double>("hltEMIsolInnerConeEB");
0064   hltEMIsolEtaSliceEB_ = conf.getParameter<double>("hltEMIsolEtaSliceEB");
0065   hltEMIsolEtMinEB_ = conf.getParameter<double>("hltEMIsolEtMinEB");
0066   hltEMIsolEMinEB_ = conf.getParameter<double>("hltEMIsolEMinEB");
0067   hltEMIsolInnerConeEE_ = conf.getParameter<double>("hltEMIsolInnerConeEE");
0068   hltEMIsolEtaSliceEE_ = conf.getParameter<double>("hltEMIsolEtaSliceEE");
0069   hltEMIsolEtMinEE_ = conf.getParameter<double>("hltEMIsolEtMinEE");
0070   hltEMIsolEMinEE_ = conf.getParameter<double>("hltEMIsolEMinEE");
0071 
0072   hltPhoTrkIsolPtMin_ = conf.getParameter<double>("hltPhoTrkIsolPtMin");
0073   hltPhoTrkIsolOuterCone_ = conf.getParameter<double>("hltPhoTrkIsolOuterCone");
0074   hltPhoTrkIsolInnerCone_ = conf.getParameter<double>("hltPhoTrkIsolInnerCone");
0075   hltPhoTrkIsolZSpan_ = conf.getParameter<double>("hltPhoTrkIsolZSpan");
0076   hltPhoTrkIsolRSpan_ = conf.getParameter<double>("hltPhoTrkIsolZSpan");
0077   hltPhoTrkIsolCountTrks_ = conf.getParameter<bool>("hltPhoTrkIsolCountTrks");
0078 
0079   hltEleTrkIsolPtMin_ = conf.getParameter<double>("hltEleTrkIsolPtMin");
0080   hltEleTrkIsolOuterCone_ = conf.getParameter<double>("hltEleTrkIsolOuterCone");
0081   hltEleTrkIsolInnerCone_ = conf.getParameter<double>("hltEleTrkIsolInnerCone");
0082   hltEleTrkIsolZSpan_ = conf.getParameter<double>("hltEleTrkIsolZSpan");
0083   hltEleTrkIsolRSpan_ = conf.getParameter<double>("hltEleTrkIsolZSpan");
0084 
0085   hltHadIsolOuterCone_ = conf.getParameter<double>("hltHadIsolOuterCone");
0086   hltHadIsolInnerCone_ = conf.getParameter<double>("hltHadIsolInnerCone");
0087   hltHadIsolEtMin_ = conf.getParameter<double>("hltHadIsolEtMin");
0088   hltHadIsolDepth_ = conf.getParameter<int>("hltHadIsolDepth");
0089 
0090   calHLTHcalIsol_ = conf.getParameter<bool>("calHLTHcalIsol");
0091   calHLTEmIsol_ = conf.getParameter<bool>("calHLTEmIsol");
0092   calHLTEleTrkIsol_ = conf.getParameter<bool>("calHLTEleTrkIsol");
0093   calHLTPhoTrkIsol_ = conf.getParameter<bool>("calHLTPhoTrkIsol");
0094 
0095   trigCutParams_ = conf.getParameter<std::vector<edm::ParameterSet> >(
0096       "triggerCuts");  //setupTriggers used to be in this function but had to be moved due to HLTConfigChanges (has to be called beginRun) so we have to save this for later.
0097 
0098   hltEleTrkIsolAlgo_ = new EgammaHLTTrackIsolation(
0099       hltEleTrkIsolPtMin_, hltEleTrkIsolOuterCone_, hltEleTrkIsolZSpan_, hltEleTrkIsolRSpan_, hltEleTrkIsolInnerCone_);
0100   hltPhoTrkIsolAlgo_ = new EgammaHLTTrackIsolation(
0101       hltPhoTrkIsolPtMin_, hltPhoTrkIsolOuterCone_, hltPhoTrkIsolZSpan_, hltPhoTrkIsolRSpan_, hltPhoTrkIsolInnerCone_);
0102 }
0103 
0104 //this code was taken out of OffHelper::setup due to HLTConfigProvider changes
0105 //it still assumes that this is called only once
0106 void OffHelper::setupTriggers(const HLTConfigProvider& hltConfig,
0107                               const std::vector<std::string>& hltFiltersUsed,
0108                               const TrigCodes& trigCodes) {
0109   hltFiltersUsed_ = hltFiltersUsed;  //expensive but only do this once and faster ways could make things less clear
0110   //now work out how many objects are requires to pass filter for it to accept
0111   hltFiltersUsedWithNrCandsCut_.clear();
0112   std::vector<int> getMRObjs = egHLT::trigTools::getMinNrObjsRequiredByFilter(hltFiltersUsed_);
0113   for (size_t filterNr = 0; filterNr < hltFiltersUsed_.size(); filterNr++) {
0114     hltFiltersUsedWithNrCandsCut_.push_back(std::make_pair(hltFiltersUsed_[filterNr], getMRObjs[filterNr]));
0115   }
0116 
0117   //now loading the cuts for every trigger into our vector which stores them
0118   //only load cuts for triggers that are in hltFiltersUsed
0119 
0120   for (auto& trigCutParam : trigCutParams_) {
0121     std::string trigName = trigCutParam.getParameter<std::string>("trigName");
0122     if (std::find(hltFiltersUsed_.begin(), hltFiltersUsed_.end(), trigName) !=
0123         hltFiltersUsed_.end()) {  //perhaps I should sort hltFiltersUsed_....
0124       trigCuts_.push_back(std::make_pair(trigCodes.getCode(trigName), OffEgSel(trigCutParam)));
0125       //   std::cout<<trigName<<std::endl<<"between"<<std::endl<<trigCutParams_[trigNr]<<std::endl<<"after"<<std::endl;
0126     }
0127   }
0128   trigCutParams_.clear();  //dont need it any more, get rid of it
0129 
0130   //to make my life difficult, the scaled l1 paths are special
0131   //and arent stored in trigger event
0132   //to I have to figure out the path, see if it passes
0133   //and then hunt down the l1 seed filter and use that to match to the pho/ele
0134   //matching on l1 seed filter is not enough as that will be passed for normal
0135   //electron triggers even if pre-scale hasnt fired
0136   l1PreScaledFilters_.clear();
0137   l1PreScaledPaths_.clear();
0138   l1PreAndSeedFilters_.clear();
0139   for (auto& filterNr : hltFiltersUsed_) {
0140     if (filterNr.find("hltPreL1") == 0) {  //l1 prescaled path
0141       l1PreScaledFilters_.push_back(filterNr);
0142     }
0143   }
0144 
0145   egHLT::trigTools::translateFiltersToPathNames(hltConfig, l1PreScaledFilters_, l1PreScaledPaths_);
0146   if (l1PreScaledPaths_.size() == l1PreScaledFilters_.size()) {
0147     for (size_t pathNr = 0; pathNr < l1PreScaledPaths_.size(); pathNr++) {
0148       std::string l1SeedFilter = egHLT::trigTools::getL1SeedFilterOfPath(hltConfig, l1PreScaledPaths_[pathNr]);
0149       //---Morse====
0150       //std::cout<<l1PreScaledFilters_[pathNr]<<"  "<<l1PreScaledPaths_[pathNr]<<"  "<<l1SeedFilter<<std::endl;
0151       //------------
0152       l1PreAndSeedFilters_.push_back(std::make_pair(l1PreScaledFilters_[pathNr], l1SeedFilter));
0153     }
0154   }
0155 }
0156 
0157 int OffHelper::makeOffEvt(const edm::Event& edmEvent,
0158                           const edm::EventSetup& setup,
0159                           egHLT::OffEvt& offEvent,
0160                           const TrigCodes& c) {
0161   offEvent.clear();
0162   int errCode = 0;  //excution stops as soon as an error is flagged
0163   if (errCode == 0)
0164     errCode = getHandles(edmEvent, setup);
0165   if (errCode == 0)
0166     errCode = fillOffEleVec(offEvent.eles());
0167   if (errCode == 0)
0168     errCode = fillOffPhoVec(offEvent.phos());
0169   if (errCode == 0)
0170     errCode = setTrigInfo(edmEvent, offEvent, c);
0171   if (errCode == 0)
0172     offEvent.setJets(recoJets_);
0173   return errCode;
0174 }
0175 
0176 int OffHelper::getHandles(const edm::Event& event, const edm::EventSetup& setup) {
0177   try {
0178     caloGeom_ = setup.getHandle(caloGeomToken_);
0179     caloTopology_ = setup.getHandle(caloTopoToken_);
0180     //ecalSeverityLevel_ = setup.getHandle(ecalSeverityToken_);
0181   } catch (cms::Exception& iException) {
0182     return errCodes::Geom;
0183   }
0184   try {
0185     magField_ = setup.getHandle(magFieldToken_);
0186   } catch (cms::Exception& iException) {
0187     return errCodes::MagField;
0188   }
0189 
0190   //get objects
0191   if (!getHandle(event, triggerSummaryToken, trigEvt_))
0192     return errCodes::TrigEvent;  //must have this, otherwise skip event
0193   if (!getHandle(event, trigResultsToken, trigResults_))
0194     return errCodes::TrigEvent;  //re using bit to minimise bug fix code changes
0195   if (!getHandle(event, electronsToken, recoEles_))
0196     return errCodes::OffEle;  //need for electrons
0197   if (!getHandle(event, photonsToken, recoPhos_))
0198     return errCodes::OffPho;  //need for photons
0199   if (!getHandle(event, caloJetsToken, recoJets_))
0200     return errCodes::OffJet;  //need for electrons and photons
0201   if (!getHandle(event, vertexToken, recoVertices_))
0202     return errCodes::OffVertex;  //need for eff vs nVertex
0203 
0204   //need for HLT isolations (rec hits also need for sigmaIPhiIPhi (ele/pho) and r9 pho)
0205   if (!getHandle(event, ecalRecHitsEBToken, ebRecHits_))
0206     return errCodes::EBRecHits;
0207   if (!getHandle(event, ecalRecHitsEEToken, eeRecHits_))
0208     return errCodes::EERecHits;
0209   if (!getHandle(event, isolTrkToken, isolTrks_))
0210     return errCodes::IsolTrks;
0211   if (!getHandle(event, hbheHitsToken, hbheHits_))
0212     return errCodes::HBHERecHits;  //I dont think we need hbhe rec-hits any more
0213   if (!getHandle(event, hfHitsToken, hfHits_))
0214     return errCodes::HFRecHits;  //I dont think we need hf rec-hits any more
0215   if (!getHandle(event, beamSpotToken, beamSpot_))
0216     return errCodes::BeamSpot;
0217   if (!getHandle(event, caloTowersToken, caloTowers_))
0218     return errCodes::CaloTowers;
0219 
0220   return 0;
0221 }
0222 
0223 //this function coverts GsfElectrons to a format which is actually useful to me
0224 int OffHelper::fillOffEleVec(std::vector<OffEle>& egHLTOffEles) {
0225   egHLTOffEles.clear();
0226   egHLTOffEles.reserve(recoEles_->size());
0227   for (auto const& gsfIter : *recoEles_) {
0228     if (!gsfIter.ecalDrivenSeed())
0229       continue;  //avoid PF electrons (this is Eg HLT validation and HLT is ecal driven)
0230 
0231     int nVertex = 0;
0232     for (auto const& nVit : *recoVertices_) {
0233       if (!nVit.isFake() && nVit.ndof() > 4 && std::fabs(nVit.z() < 24.0) &&
0234           sqrt(nVit.x() * nVit.x() + nVit.y() * nVit.y()) < 2.0) {
0235         nVertex++;
0236       }
0237     }
0238     //if(nVertex>20)std::cout<<"nVertex: "<<nVertex<<std::endl;
0239     OffEle::EventData eventData;
0240     eventData.NVertex = nVertex;
0241 
0242     OffEle::IsolData isolData;
0243     fillIsolData(gsfIter, isolData);
0244 
0245     OffEle::ClusShapeData clusShapeData;
0246     fillClusShapeData(gsfIter, clusShapeData);
0247 
0248     OffEle::HLTData hltData;
0249     fillHLTData(gsfIter, hltData);
0250 
0251     egHLTOffEles.emplace_back(gsfIter, clusShapeData, isolData, hltData, eventData);
0252 
0253     //now we would like to set the cut results
0254     OffEle& ele = egHLTOffEles.back();
0255     ele.setCutCode(eleCuts_.getCutCode(ele));
0256     ele.setLooseCutCode(eleLooseCuts_.getCutCode(ele));
0257 
0258     std::vector<std::pair<TrigCodes::TrigBitSet, int> > trigCutsCutCodes;
0259     for (auto& trigCut : trigCuts_)
0260       trigCutsCutCodes.push_back(std::make_pair(trigCut.first, trigCut.second.getCutCode(ele)));
0261     ele.setTrigCutsCutCodes(trigCutsCutCodes);
0262   }  //end loop over gsf electron collection
0263   return 0;
0264 }
0265 
0266 void OffHelper::fillIsolData(const reco::GsfElectron& ele, OffEle::IsolData& isolData) {
0267   EgammaTowerIsolation hcalIsolAlgo(
0268       hltHadIsolOuterCone_, hltHadIsolInnerCone_, hltHadIsolEtMin_, hltHadIsolDepth_, caloTowers_.product());
0269   EgammaRecHitIsolation ecalIsolAlgoEB(hltEMIsolOuterCone_,
0270                                        hltEMIsolInnerConeEB_,
0271                                        hltEMIsolEtaSliceEB_,
0272                                        hltEMIsolEtMinEB_,
0273                                        hltEMIsolEMinEB_,
0274                                        caloGeom_,
0275                                        *ebRecHits_,
0276                                        ecalSeverityLevel_.product(),
0277                                        DetId::Ecal);
0278   EgammaRecHitIsolation ecalIsolAlgoEE(hltEMIsolOuterCone_,
0279                                        hltEMIsolInnerConeEE_,
0280                                        hltEMIsolEtaSliceEE_,
0281                                        hltEMIsolEtMinEE_,
0282                                        hltEMIsolEMinEE_,
0283                                        caloGeom_,
0284                                        *eeRecHits_,
0285                                        ecalSeverityLevel_.product(),
0286                                        DetId::Ecal);
0287 
0288   isolData.ptTrks = ele.dr03TkSumPt();
0289   isolData.nrTrks = 999;  //no longer supported
0290   isolData.em = ele.dr03EcalRecHitSumEt();
0291   isolData.hadDepth1 = ele.dr03HcalTowerSumEt(1);
0292   isolData.hadDepth2 = ele.dr03HcalTowerSumEt(2);
0293 
0294   //now time to do the HLT algos
0295   if (calHLTHcalIsol_)
0296     isolData.hltHad = hcalIsolAlgo.getTowerESum(&ele);
0297   else
0298     isolData.hltHad = 0.;
0299   if (calHLTEleTrkIsol_)
0300     isolData.hltTrksEle = hltEleTrkIsolAlgo_->electronPtSum(&(*(ele.gsfTrack())), isolTrks_.product());
0301   else
0302     isolData.hltTrksEle = 0.;
0303   if (calHLTPhoTrkIsol_) {
0304     if (hltPhoTrkIsolCountTrks_)
0305       isolData.hltTrksPho = hltPhoTrkIsolAlgo_->photonTrackCount(&ele, isolTrks_.product(), false);
0306     else
0307       isolData.hltTrksPho = hltPhoTrkIsolAlgo_->photonPtSum(&ele, isolTrks_.product(), false);
0308   } else
0309     isolData.hltTrksPho = 0.;
0310   if (calHLTEmIsol_)
0311     isolData.hltEm = ecalIsolAlgoEB.getEtSum(&ele) + ecalIsolAlgoEE.getEtSum(&ele);
0312   else
0313     isolData.hltEm = 0.;
0314 }
0315 
0316 void OffHelper::fillClusShapeData(const reco::GsfElectron& ele, OffEle::ClusShapeData& clusShapeData) {
0317   clusShapeData.sigmaEtaEta = ele.sigmaEtaEta();
0318   clusShapeData.sigmaIEtaIEta = ele.sigmaIetaIeta();
0319   double e5x5 = ele.e5x5();
0320   if (e5x5 != 0.) {
0321     clusShapeData.e1x5Over5x5 = ele.e1x5() / e5x5;
0322     clusShapeData.e2x5MaxOver5x5 = ele.e2x5Max() / e5x5;
0323   } else {
0324     clusShapeData.e1x5Over5x5 = -1;
0325     clusShapeData.e2x5MaxOver5x5 = -1;
0326   }
0327 
0328   //want to calculate r9, sigmaPhiPhi and sigmaIPhiIPhi, have to do old fashioned way
0329   const reco::BasicCluster& seedClus = *(ele.superCluster()->seed());
0330   const DetId seedDetId =
0331       seedClus.hitsAndFractions()[0]
0332           .first;  //note this may not actually be the seed hit but it doesnt matter because all hits will be in the barrel OR endcap
0333   if (seedDetId.subdetId() == EcalBarrel) {
0334     const auto& stdCov =
0335         EcalClusterTools::covariances(seedClus, ebRecHits_.product(), caloTopology_.product(), caloGeom_.product());
0336     const auto& crysCov = EcalClusterTools::localCovariances(seedClus, ebRecHits_.product(), caloTopology_.product());
0337     clusShapeData.sigmaPhiPhi = sqrt(stdCov[2]);
0338     clusShapeData.sigmaIPhiIPhi = sqrt(crysCov[2]);
0339     if (ele.superCluster()->rawEnergy() != 0.) {
0340       clusShapeData.r9 = EcalClusterTools::e3x3(seedClus, ebRecHits_.product(), caloTopology_.product()) /
0341                          ele.superCluster()->rawEnergy();
0342     } else
0343       clusShapeData.r9 = -1.;
0344 
0345   } else {
0346     const auto& stdCov =
0347         EcalClusterTools::covariances(seedClus, eeRecHits_.product(), caloTopology_.product(), caloGeom_.product());
0348     const auto& crysCov = EcalClusterTools::localCovariances(seedClus, eeRecHits_.product(), caloTopology_.product());
0349     clusShapeData.sigmaPhiPhi = sqrt(stdCov[2]);
0350     clusShapeData.sigmaIPhiIPhi = sqrt(crysCov[2]);
0351     if (ele.superCluster()->rawEnergy() != 0.) {
0352       clusShapeData.r9 = EcalClusterTools::e3x3(seedClus, eeRecHits_.product(), caloTopology_.product()) /
0353                          ele.superCluster()->rawEnergy();
0354     } else
0355       clusShapeData.r9 = -1.;
0356   }
0357 }
0358 
0359 //reco approximations of hlt quantities
0360 void OffHelper::fillHLTData(const reco::GsfElectron& ele, OffEle::HLTData& hltData) {
0361   if (ele.closestCtfTrackRef().isNonnull() && ele.closestCtfTrackRef()->extra().isNonnull()) {
0362     reco::TrackRef ctfTrack = ele.closestCtfTrackRef();
0363     reco::SuperClusterRef scClus = ele.superCluster();
0364 
0365     //dEta
0366     const reco::BeamSpot::Point& bsPos = beamSpot_->position();
0367     math::XYZPoint scPosWRTVtx(scClus->x() - bsPos.x(), scClus->y() - bsPos.y(), scClus->z() - ctfTrack->vz());
0368     hltData.dEtaIn = fabs(scPosWRTVtx.eta() - ctfTrack->eta());
0369 
0370     //dPhi: lifted straight from hlt code
0371     float deltaPhi = fabs(ctfTrack->outerPosition().phi() - scClus->phi());
0372     if (deltaPhi > 6.283185308)
0373       deltaPhi -= 6.283185308;
0374     if (deltaPhi > 3.141592654)
0375       deltaPhi = 6.283185308 - deltaPhi;
0376     hltData.dPhiIn = deltaPhi;
0377 
0378     //invEInvP
0379     if (ele.ecalEnergy() != 0 && ctfTrack->p() != 0)
0380       hltData.invEInvP = 1 / ele.ecalEnergy() - 1 / ctfTrack->p();
0381     else
0382       hltData.invEInvP = 0;
0383   } else {
0384     hltData.dEtaIn = 999;
0385     hltData.dPhiIn = 999;
0386     hltData.invEInvP = 999;
0387   }
0388 
0389   //Now get HLT p4 from triggerobject
0390   //reset the position first
0391   hltData.HLTeta = 999;
0392   hltData.HLTphi = 999;
0393   hltData.HLTenergy = -999;
0394   trigTools::fillHLTposition(ele, hltData, hltFiltersUsed_, trigEvt_.product(), hltTag_);
0395   //trigTools::fillHLTposition(phos(),hltFiltersUsed_,l1PreAndSeedFilters_,evtTrigBits,trigEvt_.product(),hltTag_);
0396 }
0397 
0398 void OffHelper::fillHLTDataPho(const reco::Photon& pho, OffPho::HLTData& hltData) {
0399   //Now get HLT p4 from triggerobject
0400   //reset the position first
0401   hltData.HLTeta = 999;
0402   hltData.HLTphi = 999;
0403   hltData.HLTenergy = -999;
0404   trigTools::fillHLTposition(pho, hltData, hltFiltersUsed_, trigEvt_.product(), hltTag_);
0405   //trigTools::fillHLTposition(phos(),hltFiltersUsed_,l1PreAndSeedFilters_,evtTrigBits,trigEvt_.product(),hltTag_);
0406 }
0407 
0408 //this function coverts Photons to a format which more useful to me
0409 int OffHelper::fillOffPhoVec(std::vector<OffPho>& egHLTOffPhos) {
0410   egHLTOffPhos.clear();
0411   egHLTOffPhos.reserve(recoPhos_->size());
0412   for (auto const& phoIter : *recoPhos_) {
0413     OffPho::IsolData isolData;
0414     OffPho::ClusShapeData clusShapeData;
0415 
0416     fillIsolData(phoIter, isolData);
0417     fillClusShapeData(phoIter, clusShapeData);
0418 
0419     OffPho::HLTData hltData;
0420     fillHLTDataPho(phoIter, hltData);
0421 
0422     egHLTOffPhos.emplace_back(phoIter, clusShapeData, isolData, hltData);
0423     OffPho& pho = egHLTOffPhos.back();
0424     pho.setCutCode(phoCuts_.getCutCode(pho));
0425     pho.setLooseCutCode(phoLooseCuts_.getCutCode(pho));
0426 
0427     std::vector<std::pair<TrigCodes::TrigBitSet, int> > trigCutsCutCodes;
0428     for (auto& trigCut : trigCuts_)
0429       trigCutsCutCodes.push_back(std::make_pair(trigCut.first, trigCut.second.getCutCode(pho)));
0430     pho.setTrigCutsCutCodes(trigCutsCutCodes);
0431 
0432   }  //end loop over photon collection
0433   return 0;
0434 }
0435 
0436 void OffHelper::fillIsolData(const reco::Photon& pho, OffPho::IsolData& isolData) {
0437   EgammaTowerIsolation hcalIsolAlgo(
0438       hltHadIsolOuterCone_, hltHadIsolInnerCone_, hltHadIsolEtMin_, hltHadIsolDepth_, caloTowers_.product());
0439   EgammaRecHitIsolation ecalIsolAlgoEB(hltEMIsolOuterCone_,
0440                                        hltEMIsolInnerConeEB_,
0441                                        hltEMIsolEtaSliceEB_,
0442                                        hltEMIsolEtMinEB_,
0443                                        hltEMIsolEMinEB_,
0444                                        caloGeom_,
0445                                        *ebRecHits_,
0446                                        ecalSeverityLevel_.product(),
0447                                        DetId::Ecal);
0448   EgammaRecHitIsolation ecalIsolAlgoEE(hltEMIsolOuterCone_,
0449                                        hltEMIsolInnerConeEE_,
0450                                        hltEMIsolEtaSliceEE_,
0451                                        hltEMIsolEtMinEE_,
0452                                        hltEMIsolEMinEE_,
0453                                        caloGeom_,
0454                                        *eeRecHits_,
0455                                        ecalSeverityLevel_.product(),
0456                                        DetId::Ecal);
0457 
0458   isolData.nrTrks = pho.nTrkHollowConeDR03();
0459   isolData.ptTrks = pho.trkSumPtHollowConeDR03();
0460   isolData.em = pho.ecalRecHitSumEtConeDR03();
0461   isolData.had = pho.hcalTowerSumEtConeDR03();
0462 
0463   //now calculate hlt algos
0464   if (calHLTHcalIsol_)
0465     isolData.hltHad = hcalIsolAlgo.getTowerESum(&pho);
0466   else
0467     isolData.hltHad = 0.;
0468   if (calHLTPhoTrkIsol_) {
0469     if (hltPhoTrkIsolCountTrks_)
0470       isolData.hltTrks = hltPhoTrkIsolAlgo_->photonTrackCount(&pho, isolTrks_.product(), false);
0471     else
0472       isolData.hltTrks = hltPhoTrkIsolAlgo_->photonPtSum(&pho, isolTrks_.product(), false);
0473   } else
0474     isolData.hltTrks = 0.;
0475   if (calHLTEmIsol_)
0476     isolData.hltEm = ecalIsolAlgoEB.getEtSum(&pho) + ecalIsolAlgoEE.getEtSum(&pho);
0477   else
0478     isolData.hltEm = 0.;
0479 }
0480 
0481 void OffHelper::fillClusShapeData(const reco::Photon& pho, OffPho::ClusShapeData& clusShapeData) {
0482   clusShapeData.sigmaEtaEta = pho.sigmaEtaEta();
0483   clusShapeData.sigmaIEtaIEta = pho.sigmaIetaIeta();
0484   double e5x5 = pho.e5x5();
0485   if (e5x5 !=
0486       0.) {  //even though it is almost impossible for this to be 0., this code can never ever crash under any situation
0487     clusShapeData.e1x5Over5x5 = pho.e1x5() / e5x5;
0488     clusShapeData.e2x5MaxOver5x5 = pho.e2x5() / e5x5;
0489   } else {
0490     clusShapeData.e1x5Over5x5 = -1;
0491     clusShapeData.e2x5MaxOver5x5 = -1;
0492   }
0493   clusShapeData.r9 = pho.r9();
0494 
0495   //sigmaPhiPhi and sigmaIPhiIPhi are not in object (and nor should they be) so have to get them old fashioned way
0496   //need to figure out if its in the barrel or endcap
0497   //get the first hit of the cluster and figure out if its barrel or endcap
0498   const reco::BasicCluster& seedClus = *(pho.superCluster()->seed());
0499   const DetId seedDetId =
0500       seedClus.hitsAndFractions()[0]
0501           .first;  //note this may not actually be the seed hit but it doesnt matter because all hits will be in the barrel OR endcap (it is also incredably inefficient as it getHitsByDetId passes the vector by value not reference
0502   if (seedDetId.subdetId() == EcalBarrel) {
0503     const auto& stdCov =
0504         EcalClusterTools::covariances(seedClus, ebRecHits_.product(), caloTopology_.product(), caloGeom_.product());
0505     const auto& crysCov = EcalClusterTools::localCovariances(seedClus, ebRecHits_.product(), caloTopology_.product());
0506     clusShapeData.sigmaPhiPhi = sqrt(stdCov[2]);
0507     clusShapeData.sigmaIPhiIPhi = sqrt(crysCov[2]);
0508   } else {
0509     const auto& stdCov =
0510         EcalClusterTools::covariances(seedClus, eeRecHits_.product(), caloTopology_.product(), caloGeom_.product());
0511     const auto& crysCov = EcalClusterTools::localCovariances(seedClus, eeRecHits_.product(), caloTopology_.product());
0512     clusShapeData.sigmaPhiPhi = sqrt(stdCov[2]);
0513     clusShapeData.sigmaIPhiIPhi = sqrt(crysCov[2]);
0514   }
0515 }
0516 
0517 int OffHelper::setTrigInfo(const edm::Event& edmEvent, egHLT::OffEvt& offEvent, const TrigCodes& trigCodes) {
0518   TrigCodes::TrigBitSet evtTrigBits =
0519       trigTools::getFiltersPassed(hltFiltersUsedWithNrCandsCut_, trigEvt_.product(), hltTag_, trigCodes);
0520   //the l1 prescale paths dont have a filter with I can figure out if it passed or failed with so have to use TriggerResults
0521   if (l1PreScaledPaths_.size() ==
0522       l1PreScaledFilters_.size()) {  //check to ensure both vectors have same number of events incase of screw ups
0523     const edm::TriggerNames& triggerNames = edmEvent.triggerNames(*trigResults_);
0524     for (size_t pathNr = 0; pathNr < l1PreScaledPaths_.size();
0525          pathNr++) {  //now we have to check the prescaled l1 trigger paths
0526       unsigned int pathIndex = triggerNames.triggerIndex(l1PreScaledPaths_[pathNr]);
0527       if (pathIndex < trigResults_->size() && trigResults_->accept(pathIndex)) {
0528         evtTrigBits |= trigCodes.getCode(l1PreScaledFilters_[pathNr]);
0529       }
0530     }
0531   }
0532 
0533   offEvent.setEvtTrigBits(evtTrigBits);
0534 
0535   trigTools::setFiltersObjPasses(
0536       offEvent.eles(), hltFiltersUsed_, l1PreAndSeedFilters_, evtTrigBits, trigCodes, trigEvt_.product(), hltTag_);
0537   trigTools::setFiltersObjPasses(
0538       offEvent.phos(), hltFiltersUsed_, l1PreAndSeedFilters_, evtTrigBits, trigCodes, trigEvt_.product(), hltTag_);
0539   return 0;
0540 }
[ 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!