Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-10-08 05:11:33 

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     trigCutsCutCodes.reserve(trigCuts_.size());
0260     for (auto& trigCut : trigCuts_)
0261       trigCutsCutCodes.push_back(std::make_pair(trigCut.first, trigCut.second.getCutCode(ele)));
0262     ele.setTrigCutsCutCodes(trigCutsCutCodes);
0263   }  //end loop over gsf electron collection
0264   return 0;
0265 }
0266 
0267 void OffHelper::fillIsolData(const reco::GsfElectron& ele, OffEle::IsolData& isolData) {
0268   EgammaTowerIsolation hcalIsolAlgo(
0269       hltHadIsolOuterCone_, hltHadIsolInnerCone_, hltHadIsolEtMin_, hltHadIsolDepth_, caloTowers_.product());
0270   EgammaRecHitIsolation ecalIsolAlgoEB(hltEMIsolOuterCone_,
0271                                        hltEMIsolInnerConeEB_,
0272                                        hltEMIsolEtaSliceEB_,
0273                                        hltEMIsolEtMinEB_,
0274                                        hltEMIsolEMinEB_,
0275                                        caloGeom_,
0276                                        *ebRecHits_,
0277                                        ecalSeverityLevel_.product(),
0278                                        DetId::Ecal);
0279   EgammaRecHitIsolation ecalIsolAlgoEE(hltEMIsolOuterCone_,
0280                                        hltEMIsolInnerConeEE_,
0281                                        hltEMIsolEtaSliceEE_,
0282                                        hltEMIsolEtMinEE_,
0283                                        hltEMIsolEMinEE_,
0284                                        caloGeom_,
0285                                        *eeRecHits_,
0286                                        ecalSeverityLevel_.product(),
0287                                        DetId::Ecal);
0288 
0289   isolData.ptTrks = ele.dr03TkSumPt();
0290   isolData.nrTrks = 999;  //no longer supported
0291   isolData.em = ele.dr03EcalRecHitSumEt();
0292   isolData.hadDepth1 = ele.dr03HcalTowerSumEt(1);
0293   isolData.hadDepth2 = ele.dr03HcalTowerSumEt(2);
0294 
0295   //now time to do the HLT algos
0296   if (calHLTHcalIsol_)
0297     isolData.hltHad = hcalIsolAlgo.getTowerESum(&ele);
0298   else
0299     isolData.hltHad = 0.;
0300   if (calHLTEleTrkIsol_)
0301     isolData.hltTrksEle = hltEleTrkIsolAlgo_->electronPtSum(&(*(ele.gsfTrack())), isolTrks_.product());
0302   else
0303     isolData.hltTrksEle = 0.;
0304   if (calHLTPhoTrkIsol_) {
0305     if (hltPhoTrkIsolCountTrks_)
0306       isolData.hltTrksPho = hltPhoTrkIsolAlgo_->photonTrackCount(&ele, isolTrks_.product(), false);
0307     else
0308       isolData.hltTrksPho = hltPhoTrkIsolAlgo_->photonPtSum(&ele, isolTrks_.product(), false);
0309   } else
0310     isolData.hltTrksPho = 0.;
0311   if (calHLTEmIsol_)
0312     isolData.hltEm = ecalIsolAlgoEB.getEtSum(&ele, thresholds) + ecalIsolAlgoEE.getEtSum(&ele, thresholds);
0313   else
0314     isolData.hltEm = 0.;
0315 }
0316 
0317 void OffHelper::fillClusShapeData(const reco::GsfElectron& ele, OffEle::ClusShapeData& clusShapeData) {
0318   clusShapeData.sigmaEtaEta = ele.sigmaEtaEta();
0319   clusShapeData.sigmaIEtaIEta = ele.sigmaIetaIeta();
0320   double e5x5 = ele.e5x5();
0321   if (e5x5 != 0.) {
0322     clusShapeData.e1x5Over5x5 = ele.e1x5() / e5x5;
0323     clusShapeData.e2x5MaxOver5x5 = ele.e2x5Max() / e5x5;
0324   } else {
0325     clusShapeData.e1x5Over5x5 = -1;
0326     clusShapeData.e2x5MaxOver5x5 = -1;
0327   }
0328 
0329   //want to calculate r9, sigmaPhiPhi and sigmaIPhiIPhi, have to do old fashioned way
0330   const reco::BasicCluster& seedClus = *(ele.superCluster()->seed());
0331   const DetId seedDetId =
0332       seedClus.hitsAndFractions()[0]
0333           .first;  //note this may not actually be the seed hit but it doesnt matter because all hits will be in the barrel OR endcap
0334   if (seedDetId.subdetId() == EcalBarrel) {
0335     const auto& stdCov =
0336         EcalClusterTools::covariances(seedClus, ebRecHits_.product(), caloTopology_.product(), caloGeom_.product());
0337     const auto& crysCov = EcalClusterTools::localCovariances(seedClus, ebRecHits_.product(), caloTopology_.product());
0338     clusShapeData.sigmaPhiPhi = sqrt(stdCov[2]);
0339     clusShapeData.sigmaIPhiIPhi = sqrt(crysCov[2]);
0340     if (ele.superCluster()->rawEnergy() != 0.) {
0341       clusShapeData.r9 = EcalClusterTools::e3x3(seedClus, ebRecHits_.product(), caloTopology_.product()) /
0342                          ele.superCluster()->rawEnergy();
0343     } else
0344       clusShapeData.r9 = -1.;
0345 
0346   } else {
0347     const auto& stdCov =
0348         EcalClusterTools::covariances(seedClus, eeRecHits_.product(), caloTopology_.product(), caloGeom_.product());
0349     const auto& crysCov = EcalClusterTools::localCovariances(seedClus, eeRecHits_.product(), caloTopology_.product());
0350     clusShapeData.sigmaPhiPhi = sqrt(stdCov[2]);
0351     clusShapeData.sigmaIPhiIPhi = sqrt(crysCov[2]);
0352     if (ele.superCluster()->rawEnergy() != 0.) {
0353       clusShapeData.r9 = EcalClusterTools::e3x3(seedClus, eeRecHits_.product(), caloTopology_.product()) /
0354                          ele.superCluster()->rawEnergy();
0355     } else
0356       clusShapeData.r9 = -1.;
0357   }
0358 }
0359 
0360 //reco approximations of hlt quantities
0361 void OffHelper::fillHLTData(const reco::GsfElectron& ele, OffEle::HLTData& hltData) {
0362   if (ele.closestCtfTrackRef().isNonnull() && ele.closestCtfTrackRef()->extra().isNonnull()) {
0363     reco::TrackRef ctfTrack = ele.closestCtfTrackRef();
0364     reco::SuperClusterRef scClus = ele.superCluster();
0365 
0366     //dEta
0367     const reco::BeamSpot::Point& bsPos = beamSpot_->position();
0368     math::XYZPoint scPosWRTVtx(scClus->x() - bsPos.x(), scClus->y() - bsPos.y(), scClus->z() - ctfTrack->vz());
0369     hltData.dEtaIn = fabs(scPosWRTVtx.eta() - ctfTrack->eta());
0370 
0371     //dPhi: lifted straight from hlt code
0372     float deltaPhi = fabs(ctfTrack->outerPosition().phi() - scClus->phi());
0373     if (deltaPhi > 6.283185308)
0374       deltaPhi -= 6.283185308;
0375     if (deltaPhi > 3.141592654)
0376       deltaPhi = 6.283185308 - deltaPhi;
0377     hltData.dPhiIn = deltaPhi;
0378 
0379     //invEInvP
0380     if (ele.ecalEnergy() != 0 && ctfTrack->p() != 0)
0381       hltData.invEInvP = 1 / ele.ecalEnergy() - 1 / ctfTrack->p();
0382     else
0383       hltData.invEInvP = 0;
0384   } else {
0385     hltData.dEtaIn = 999;
0386     hltData.dPhiIn = 999;
0387     hltData.invEInvP = 999;
0388   }
0389 
0390   //Now get HLT p4 from triggerobject
0391   //reset the position first
0392   hltData.HLTeta = 999;
0393   hltData.HLTphi = 999;
0394   hltData.HLTenergy = -999;
0395   trigTools::fillHLTposition(ele, hltData, hltFiltersUsed_, trigEvt_.product(), hltTag_);
0396   //trigTools::fillHLTposition(phos(),hltFiltersUsed_,l1PreAndSeedFilters_,evtTrigBits,trigEvt_.product(),hltTag_);
0397 }
0398 
0399 void OffHelper::fillHLTDataPho(const reco::Photon& pho, OffPho::HLTData& hltData) {
0400   //Now get HLT p4 from triggerobject
0401   //reset the position first
0402   hltData.HLTeta = 999;
0403   hltData.HLTphi = 999;
0404   hltData.HLTenergy = -999;
0405   trigTools::fillHLTposition(pho, hltData, hltFiltersUsed_, trigEvt_.product(), hltTag_);
0406   //trigTools::fillHLTposition(phos(),hltFiltersUsed_,l1PreAndSeedFilters_,evtTrigBits,trigEvt_.product(),hltTag_);
0407 }
0408 
0409 //this function coverts Photons to a format which more useful to me
0410 int OffHelper::fillOffPhoVec(std::vector<OffPho>& egHLTOffPhos) {
0411   egHLTOffPhos.clear();
0412   egHLTOffPhos.reserve(recoPhos_->size());
0413   for (auto const& phoIter : *recoPhos_) {
0414     OffPho::IsolData isolData;
0415     OffPho::ClusShapeData clusShapeData;
0416 
0417     fillIsolData(phoIter, isolData);
0418     fillClusShapeData(phoIter, clusShapeData);
0419 
0420     OffPho::HLTData hltData;
0421     fillHLTDataPho(phoIter, hltData);
0422 
0423     egHLTOffPhos.emplace_back(phoIter, clusShapeData, isolData, hltData);
0424     OffPho& pho = egHLTOffPhos.back();
0425     pho.setCutCode(phoCuts_.getCutCode(pho));
0426     pho.setLooseCutCode(phoLooseCuts_.getCutCode(pho));
0427 
0428     std::vector<std::pair<TrigCodes::TrigBitSet, int> > trigCutsCutCodes;
0429     trigCutsCutCodes.reserve(trigCuts_.size());
0430     for (auto& trigCut : trigCuts_)
0431       trigCutsCutCodes.push_back(std::make_pair(trigCut.first, trigCut.second.getCutCode(pho)));
0432     pho.setTrigCutsCutCodes(trigCutsCutCodes);
0433 
0434   }  //end loop over photon collection
0435   return 0;
0436 }
0437 
0438 void OffHelper::fillIsolData(const reco::Photon& pho, OffPho::IsolData& isolData) {
0439   EgammaTowerIsolation hcalIsolAlgo(
0440       hltHadIsolOuterCone_, hltHadIsolInnerCone_, hltHadIsolEtMin_, hltHadIsolDepth_, caloTowers_.product());
0441   EgammaRecHitIsolation ecalIsolAlgoEB(hltEMIsolOuterCone_,
0442                                        hltEMIsolInnerConeEB_,
0443                                        hltEMIsolEtaSliceEB_,
0444                                        hltEMIsolEtMinEB_,
0445                                        hltEMIsolEMinEB_,
0446                                        caloGeom_,
0447                                        *ebRecHits_,
0448                                        ecalSeverityLevel_.product(),
0449                                        DetId::Ecal);
0450   EgammaRecHitIsolation ecalIsolAlgoEE(hltEMIsolOuterCone_,
0451                                        hltEMIsolInnerConeEE_,
0452                                        hltEMIsolEtaSliceEE_,
0453                                        hltEMIsolEtMinEE_,
0454                                        hltEMIsolEMinEE_,
0455                                        caloGeom_,
0456                                        *eeRecHits_,
0457                                        ecalSeverityLevel_.product(),
0458                                        DetId::Ecal);
0459 
0460   isolData.nrTrks = pho.nTrkHollowConeDR03();
0461   isolData.ptTrks = pho.trkSumPtHollowConeDR03();
0462   isolData.em = pho.ecalRecHitSumEtConeDR03();
0463   isolData.had = pho.hcalTowerSumEtConeDR03();
0464 
0465   //now calculate hlt algos
0466   if (calHLTHcalIsol_)
0467     isolData.hltHad = hcalIsolAlgo.getTowerESum(&pho);
0468   else
0469     isolData.hltHad = 0.;
0470   if (calHLTPhoTrkIsol_) {
0471     if (hltPhoTrkIsolCountTrks_)
0472       isolData.hltTrks = hltPhoTrkIsolAlgo_->photonTrackCount(&pho, isolTrks_.product(), false);
0473     else
0474       isolData.hltTrks = hltPhoTrkIsolAlgo_->photonPtSum(&pho, isolTrks_.product(), false);
0475   } else
0476     isolData.hltTrks = 0.;
0477   if (calHLTEmIsol_)
0478     isolData.hltEm = ecalIsolAlgoEB.getEtSum(&pho, thresholds) + ecalIsolAlgoEE.getEtSum(&pho, thresholds);
0479   else
0480     isolData.hltEm = 0.;
0481 }
0482 
0483 void OffHelper::fillClusShapeData(const reco::Photon& pho, OffPho::ClusShapeData& clusShapeData) {
0484   clusShapeData.sigmaEtaEta = pho.sigmaEtaEta();
0485   clusShapeData.sigmaIEtaIEta = pho.sigmaIetaIeta();
0486   double e5x5 = pho.e5x5();
0487   if (e5x5 !=
0488       0.) {  //even though it is almost impossible for this to be 0., this code can never ever crash under any situation
0489     clusShapeData.e1x5Over5x5 = pho.e1x5() / e5x5;
0490     clusShapeData.e2x5MaxOver5x5 = pho.e2x5() / e5x5;
0491   } else {
0492     clusShapeData.e1x5Over5x5 = -1;
0493     clusShapeData.e2x5MaxOver5x5 = -1;
0494   }
0495   clusShapeData.r9 = pho.r9();
0496 
0497   //sigmaPhiPhi and sigmaIPhiIPhi are not in object (and nor should they be) so have to get them old fashioned way
0498   //need to figure out if its in the barrel or endcap
0499   //get the first hit of the cluster and figure out if its barrel or endcap
0500   const reco::BasicCluster& seedClus = *(pho.superCluster()->seed());
0501   const DetId seedDetId =
0502       seedClus.hitsAndFractions()[0]
0503           .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
0504   if (seedDetId.subdetId() == EcalBarrel) {
0505     const auto& stdCov =
0506         EcalClusterTools::covariances(seedClus, ebRecHits_.product(), caloTopology_.product(), caloGeom_.product());
0507     const auto& crysCov = EcalClusterTools::localCovariances(seedClus, ebRecHits_.product(), caloTopology_.product());
0508     clusShapeData.sigmaPhiPhi = sqrt(stdCov[2]);
0509     clusShapeData.sigmaIPhiIPhi = sqrt(crysCov[2]);
0510   } else {
0511     const auto& stdCov =
0512         EcalClusterTools::covariances(seedClus, eeRecHits_.product(), caloTopology_.product(), caloGeom_.product());
0513     const auto& crysCov = EcalClusterTools::localCovariances(seedClus, eeRecHits_.product(), caloTopology_.product());
0514     clusShapeData.sigmaPhiPhi = sqrt(stdCov[2]);
0515     clusShapeData.sigmaIPhiIPhi = sqrt(crysCov[2]);
0516   }
0517 }
0518 
0519 int OffHelper::setTrigInfo(const edm::Event& edmEvent, egHLT::OffEvt& offEvent, const TrigCodes& trigCodes) {
0520   TrigCodes::TrigBitSet evtTrigBits =
0521       trigTools::getFiltersPassed(hltFiltersUsedWithNrCandsCut_, trigEvt_.product(), hltTag_, trigCodes);
0522   //the l1 prescale paths dont have a filter with I can figure out if it passed or failed with so have to use TriggerResults
0523   if (l1PreScaledPaths_.size() ==
0524       l1PreScaledFilters_.size()) {  //check to ensure both vectors have same number of events incase of screw ups
0525     const edm::TriggerNames& triggerNames = edmEvent.triggerNames(*trigResults_);
0526     for (size_t pathNr = 0; pathNr < l1PreScaledPaths_.size();
0527          pathNr++) {  //now we have to check the prescaled l1 trigger paths
0528       unsigned int pathIndex = triggerNames.triggerIndex(l1PreScaledPaths_[pathNr]);
0529       if (pathIndex < trigResults_->size() && trigResults_->accept(pathIndex)) {
0530         evtTrigBits |= trigCodes.getCode(l1PreScaledFilters_[pathNr]);
0531       }
0532     }
0533   }
0534 
0535   offEvent.setEvtTrigBits(evtTrigBits);
0536 
0537   trigTools::setFiltersObjPasses(
0538       offEvent.eles(), hltFiltersUsed_, l1PreAndSeedFilters_, evtTrigBits, trigCodes, trigEvt_.product(), hltTag_);
0539   trigTools::setFiltersObjPasses(
0540       offEvent.phos(), hltFiltersUsed_, l1PreAndSeedFilters_, evtTrigBits, trigCodes, trigEvt_.product(), hltTag_);
0541   return 0;
0542 }
[ 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!