Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:20:37 

0001 #include "CaloExtractorByAssociator.h"
0002 
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 
0007 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0008 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0009 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0010 
0011 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0012 
0013 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0014 
0015 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0016 #include "TrackingTools/TrackAssociator/interface/TrackAssociatorParameters.h"
0017 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0018 
0019 #include "DataFormats/Math/interface/deltaR.h"
0020 
0021 using namespace edm;
0022 using namespace std;
0023 using namespace reco;
0024 using namespace muonisolation;
0025 using reco::isodeposit::Direction;
0026 
0027 namespace {
0028   constexpr double dRMax_CandDep = 1.0;  //pick up candidate own deposits up to this dR if theDR_Max is smaller
0029 }
0030 
0031 CaloExtractorByAssociator::CaloExtractorByAssociator(const ParameterSet& par, edm::ConsumesCollector&& iC)
0032     : theUseRecHitsFlag(par.getParameter<bool>("UseRecHitsFlag")),
0033       theDepositLabel(par.getUntrackedParameter<string>("DepositLabel")),
0034       theDepositInstanceLabels(par.getParameter<std::vector<std::string> >("DepositInstanceLabels")),
0035       thePropagatorName(par.getParameter<std::string>("PropagatorName")),
0036       theThreshold_E(par.getParameter<double>("Threshold_E")),
0037       theThreshold_H(par.getParameter<double>("Threshold_H")),
0038       theThreshold_HO(par.getParameter<double>("Threshold_HO")),
0039       theDR_Veto_E(par.getParameter<double>("DR_Veto_E")),
0040       theDR_Veto_H(par.getParameter<double>("DR_Veto_H")),
0041       theDR_Veto_HO(par.getParameter<double>("DR_Veto_HO")),
0042       theCenterConeOnCalIntersection(par.getParameter<bool>("CenterConeOnCalIntersection")),
0043       theDR_Max(par.getParameter<double>("DR_Max")),
0044       theNoise_EB(par.getParameter<double>("Noise_EB")),
0045       theNoise_EE(par.getParameter<double>("Noise_EE")),
0046       theNoise_HB(par.getParameter<double>("Noise_HB")),
0047       theNoise_HE(par.getParameter<double>("Noise_HE")),
0048       theNoise_HO(par.getParameter<double>("Noise_HO")),
0049       theNoiseTow_EB(par.getParameter<double>("NoiseTow_EB")),
0050       theNoiseTow_EE(par.getParameter<double>("NoiseTow_EE")),
0051       theService(nullptr),
0052       theAssociator(nullptr),
0053       bFieldToken_(iC.esConsumes()),
0054       thePrintTimeReport(par.getUntrackedParameter<bool>("PrintTimeReport")) {
0055   ParameterSet serviceParameters = par.getParameter<ParameterSet>("ServiceParameters");
0056   theService = new MuonServiceProxy(serviceParameters, edm::ConsumesCollector(iC));
0057 
0058   //theAssociatorParameters = new TrackAssociatorParameters(par.getParameter<edm::ParameterSet>("TrackAssociatorParameters"), iC);
0059   theAssociatorParameters = new TrackAssociatorParameters();
0060   theAssociatorParameters->loadParameters(par.getParameter<edm::ParameterSet>("TrackAssociatorParameters"), iC);
0061   theAssociator = new TrackDetectorAssociator();
0062 
0063   if (theUseRecHitsFlag) {
0064     caloGeomToken_ = iC.esConsumes();
0065   }
0066 }
0067 
0068 CaloExtractorByAssociator::~CaloExtractorByAssociator() {
0069   if (theAssociatorParameters)
0070     delete theAssociatorParameters;
0071   if (theService)
0072     delete theService;
0073   if (theAssociator)
0074     delete theAssociator;
0075 }
0076 
0077 void CaloExtractorByAssociator::fillVetos(const edm::Event& event,
0078                                           const edm::EventSetup& eventSetup,
0079                                           const TrackCollection& muons) {
0080   //   LogWarning("CaloExtractorByAssociator")
0081   //     <<"fillVetos does nothing now: IsoDeposit provides enough functionality\n"
0082   //     <<"to remove a deposit at/around given (eta, phi)";
0083 }
0084 
0085 IsoDeposit CaloExtractorByAssociator::deposit(const Event& event,
0086                                               const EventSetup& eventSetup,
0087                                               const Track& muon) const {
0088   IsoDeposit::Direction muonDir(muon.eta(), muon.phi());
0089   IsoDeposit dep(muonDir);
0090 
0091   //   LogWarning("CaloExtractorByAssociator")
0092   //     <<"single deposit is not an option here\n"
0093   //     <<"use ::deposits --> extract all and reweight as necessary";
0094 
0095   return dep;
0096 }
0097 
0098 //! Make separate deposits: for ECAL, HCAL, HO
0099 std::vector<IsoDeposit> CaloExtractorByAssociator::deposits(const Event& event,
0100                                                             const EventSetup& eventSetup,
0101                                                             const Track& muon) const {
0102   theService->update(eventSetup);
0103   theAssociator->setPropagator(&*(theService->propagator(thePropagatorName)));
0104 
0105   //! check configuration consistency
0106   //! could've been made at construction stage (fix later?)
0107   if (theDepositInstanceLabels.size() != 3) {
0108     LogError("MuonIsolation") << "Configuration is inconsistent: Need 3 deposit instance labels";
0109   }
0110   if (!(theDepositInstanceLabels[0].compare(0, 1, std::string("e")) == 0) ||
0111       !(theDepositInstanceLabels[1].compare(0, 1, std::string("h")) == 0) ||
0112       !(theDepositInstanceLabels[2].compare(0, 2, std::string("ho")) == 0)) {
0113     LogWarning("MuonIsolation")
0114         << "Deposit instance labels do not look like  (e*, h*, ho*):"
0115         << "proceed at your own risk. The extractor interprets lab0=from ecal; lab1=from hcal; lab2=from ho";
0116   }
0117 
0118   typedef IsoDeposit::Veto Veto;
0119   //! this should be (eventually) set to the eta-phi of the crossing point of
0120   //! a straight line tangent to a muon at IP and the calorimeter
0121   IsoDeposit::Direction muonDir(muon.eta(), muon.phi());
0122 
0123   IsoDeposit depEcal(muonDir);
0124   IsoDeposit depHcal(muonDir);
0125   IsoDeposit depHOcal(muonDir);
0126 
0127   auto const& bField = eventSetup.getData(bFieldToken_);
0128 
0129   reco::TransientTrack tMuon(muon, &bField);
0130   FreeTrajectoryState iFTS = tMuon.initialFreeState();
0131   TrackDetMatchInfo mInfo = theAssociator->associate(event, eventSetup, iFTS, *theAssociatorParameters);
0132 
0133   //! each deposit type veto is at the point of intersect with that detector
0134   depEcal.setVeto(
0135       Veto(reco::isodeposit::Direction(mInfo.trkGlobPosAtEcal.eta(), mInfo.trkGlobPosAtEcal.phi()), theDR_Veto_E));
0136   depHcal.setVeto(
0137       Veto(reco::isodeposit::Direction(mInfo.trkGlobPosAtHcal.eta(), mInfo.trkGlobPosAtHcal.phi()), theDR_Veto_H));
0138   depHOcal.setVeto(
0139       Veto(reco::isodeposit::Direction(mInfo.trkGlobPosAtHO.eta(), mInfo.trkGlobPosAtHO.phi()), theDR_Veto_HO));
0140 
0141   if (theCenterConeOnCalIntersection) {
0142     reco::isodeposit::Direction dirTmp = depEcal.veto().vetoDir;
0143     double dRtmp = depEcal.veto().dR;
0144     depEcal = IsoDeposit(dirTmp);
0145     depEcal.setVeto(Veto(dirTmp, dRtmp));
0146 
0147     dirTmp = depHcal.veto().vetoDir;
0148     dRtmp = depHcal.veto().dR;
0149     depHcal = IsoDeposit(dirTmp);
0150     depHcal.setVeto(Veto(dirTmp, dRtmp));
0151 
0152     dirTmp = depHOcal.veto().vetoDir;
0153     dRtmp = depHOcal.veto().dR;
0154     depHOcal = IsoDeposit(dirTmp);
0155     depHOcal.setVeto(Veto(dirTmp, dRtmp));
0156   }
0157 
0158   if (theUseRecHitsFlag) {
0159     //! do things based on rec-hits here
0160     //! too much copy-pasting now (refactor later?)
0161     auto const& caloGeom = eventSetup.getData(caloGeomToken_);
0162 
0163     //Ecal
0164     std::vector<const EcalRecHit*>::const_iterator eHitCI = mInfo.ecalRecHits.begin();
0165     for (; eHitCI != mInfo.ecalRecHits.end(); ++eHitCI) {
0166       const EcalRecHit* eHitCPtr = *eHitCI;
0167       GlobalPoint eHitPos = caloGeom.getPosition(eHitCPtr->detid());
0168       double deltar0 = reco::deltaR(muon, eHitPos);
0169       double cosTheta = 1. / cosh(eHitPos.eta());
0170       double energy = eHitCPtr->energy();
0171       double et = energy * cosTheta;
0172       if (deltar0 > std::max(dRMax_CandDep, theDR_Max) ||
0173           !(et > theThreshold_E && energy > 3 * noiseRecHit(eHitCPtr->detid())))
0174         continue;
0175 
0176       bool vetoHit = false;
0177       double deltar = reco::deltaR(mInfo.trkGlobPosAtEcal, eHitPos);
0178       //! first check if the hit is inside the veto cone by dR-alone
0179       if (deltar < theDR_Veto_E) {
0180         LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto ECAL hit: Calo deltaR= " << deltar;
0181         LogDebug("RecoMuon|CaloExtractorByAssociator")
0182             << " >>> Calo eta phi ethcal: " << eHitPos.eta() << " " << eHitPos.phi() << " " << et;
0183         vetoHit = true;
0184       }
0185       //! and now pitch those in the crossed list
0186       if (!vetoHit) {
0187         for (unsigned int iH = 0; iH < mInfo.crossedEcalIds.size() && !vetoHit; ++iH) {
0188           if (mInfo.crossedEcalIds[iH].rawId() == eHitCPtr->detid().rawId())
0189             vetoHit = true;
0190         }
0191       }
0192 
0193       //check theDR_Max only here to keep vetoHits being added to the veto energy
0194       if (deltar0 > theDR_Max && !vetoHit)
0195         continue;
0196 
0197       if (vetoHit) {
0198         depEcal.addCandEnergy(et);
0199       } else {
0200         depEcal.addDeposit(reco::isodeposit::Direction(eHitPos.eta(), eHitPos.phi()), et);
0201       }
0202     }
0203 
0204     //Hcal
0205     std::vector<const HBHERecHit*>::const_iterator hHitCI = mInfo.hcalRecHits.begin();
0206     for (; hHitCI != mInfo.hcalRecHits.end(); ++hHitCI) {
0207       const HBHERecHit* hHitCPtr = *hHitCI;
0208       GlobalPoint hHitPos = caloGeom.getPosition(hHitCPtr->detid());
0209       double deltar0 = reco::deltaR(muon, hHitPos);
0210       double cosTheta = 1. / cosh(hHitPos.eta());
0211       double energy = hHitCPtr->energy();
0212       double et = energy * cosTheta;
0213       if (deltar0 > std::max(dRMax_CandDep, theDR_Max) ||
0214           !(et > theThreshold_H && energy > 3 * noiseRecHit(hHitCPtr->detid())))
0215         continue;
0216 
0217       bool vetoHit = false;
0218       double deltar = reco::deltaR(mInfo.trkGlobPosAtHcal, hHitPos);
0219       //! first check if the hit is inside the veto cone by dR-alone
0220       if (deltar < theDR_Veto_H) {
0221         LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto HBHE hit: Calo deltaR= " << deltar;
0222         LogDebug("RecoMuon|CaloExtractorByAssociator")
0223             << " >>> Calo eta phi ethcal: " << hHitPos.eta() << " " << hHitPos.phi() << " " << et;
0224         vetoHit = true;
0225       }
0226       //! and now pitch those in the crossed list
0227       if (!vetoHit) {
0228         for (unsigned int iH = 0; iH < mInfo.crossedHcalIds.size() && !vetoHit; ++iH) {
0229           if (mInfo.crossedHcalIds[iH].rawId() == hHitCPtr->detid().rawId())
0230             vetoHit = true;
0231         }
0232       }
0233 
0234       //check theDR_Max only here to keep vetoHits being added to the veto energy
0235       if (deltar0 > theDR_Max && !vetoHit)
0236         continue;
0237 
0238       if (vetoHit) {
0239         depHcal.addCandEnergy(et);
0240       } else {
0241         depHcal.addDeposit(reco::isodeposit::Direction(hHitPos.eta(), hHitPos.phi()), et);
0242       }
0243     }
0244 
0245     //HOcal
0246     std::vector<const HORecHit*>::const_iterator hoHitCI = mInfo.hoRecHits.begin();
0247     for (; hoHitCI != mInfo.hoRecHits.end(); ++hoHitCI) {
0248       const HORecHit* hoHitCPtr = *hoHitCI;
0249       GlobalPoint hoHitPos = caloGeom.getPosition(hoHitCPtr->detid());
0250       double deltar0 = reco::deltaR(muon, hoHitPos);
0251       double cosTheta = 1. / cosh(hoHitPos.eta());
0252       double energy = hoHitCPtr->energy();
0253       double et = energy * cosTheta;
0254       if (deltar0 > std::max(dRMax_CandDep, theDR_Max) ||
0255           !(et > theThreshold_HO && energy > 3 * noiseRecHit(hoHitCPtr->detid())))
0256         continue;
0257 
0258       bool vetoHit = false;
0259       double deltar = reco::deltaR(mInfo.trkGlobPosAtHO, hoHitPos);
0260       //! first check if the hit is inside the veto cone by dR-alone
0261       if (deltar < theDR_Veto_HO) {
0262         LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto HO hit: Calo deltaR= " << deltar;
0263         LogDebug("RecoMuon|CaloExtractorByAssociator")
0264             << " >>> Calo eta phi ethcal: " << hoHitPos.eta() << " " << hoHitPos.phi() << " " << et;
0265         vetoHit = true;
0266       }
0267       //! and now pitch those in the crossed list
0268       if (!vetoHit) {
0269         for (unsigned int iH = 0; iH < mInfo.crossedHOIds.size() && !vetoHit; ++iH) {
0270           if (mInfo.crossedHOIds[iH].rawId() == hoHitCPtr->detid().rawId())
0271             vetoHit = true;
0272         }
0273       }
0274 
0275       //check theDR_Max only here to keep vetoHits being added to the veto energy
0276       if (deltar0 > theDR_Max && !vetoHit)
0277         continue;
0278 
0279       if (vetoHit) {
0280         depHOcal.addCandEnergy(et);
0281       } else {
0282         depHOcal.addDeposit(reco::isodeposit::Direction(hoHitPos.eta(), hoHitPos.phi()), et);
0283       }
0284     }
0285 
0286   } else {
0287     //! use calo towers
0288     std::vector<const CaloTower*>::const_iterator calCI = mInfo.towers.begin();
0289     for (; calCI != mInfo.towers.end(); ++calCI) {
0290       const CaloTower* calCPtr = *calCI;
0291       double deltar0 = reco::deltaR(muon, *calCPtr);
0292       if (deltar0 > std::max(dRMax_CandDep, theDR_Max))
0293         continue;
0294 
0295       //even more copy-pasting .. need to refactor
0296       double etecal = calCPtr->emEt();
0297       double eecal = calCPtr->emEnergy();
0298       bool doEcal = etecal > theThreshold_E && eecal > 3 * noiseEcal(*calCPtr);
0299       double ethcal = calCPtr->hadEt();
0300       double ehcal = calCPtr->hadEnergy();
0301       bool doHcal = ethcal > theThreshold_H && ehcal > 3 * noiseHcal(*calCPtr);
0302       double ethocal = calCPtr->outerEt();
0303       double ehocal = calCPtr->outerEnergy();
0304       bool doHOcal = ethocal > theThreshold_HO && ehocal > 3 * noiseHOcal(*calCPtr);
0305       if ((!doEcal) && (!doHcal) && (!doHcal))
0306         continue;
0307 
0308       bool vetoTowerEcal = false;
0309       double deltarEcal = reco::deltaR(mInfo.trkGlobPosAtEcal, *calCPtr);
0310       //! first check if the tower is inside the veto cone by dR-alone
0311       if (deltarEcal < theDR_Veto_E) {
0312         LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto ecal tower: Calo deltaR= " << deltarEcal;
0313         LogDebug("RecoMuon|CaloExtractorByAssociator")
0314             << " >>> Calo eta phi ethcal: " << calCPtr->eta() << " " << calCPtr->phi() << " " << ethcal;
0315         vetoTowerEcal = true;
0316       }
0317       bool vetoTowerHcal = false;
0318       double deltarHcal = reco::deltaR(mInfo.trkGlobPosAtHcal, *calCPtr);
0319       //! first check if the tower is inside the veto cone by dR-alone
0320       if (deltarHcal < theDR_Veto_H) {
0321         LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto hcal tower: Calo deltaR= " << deltarHcal;
0322         LogDebug("RecoMuon|CaloExtractorByAssociator")
0323             << " >>> Calo eta phi ethcal: " << calCPtr->eta() << " " << calCPtr->phi() << " " << ethcal;
0324         vetoTowerHcal = true;
0325       }
0326       bool vetoTowerHOCal = false;
0327       double deltarHOcal = reco::deltaR(mInfo.trkGlobPosAtHO, *calCPtr);
0328       //! first check if the tower is inside the veto cone by dR-alone
0329       if (deltarHOcal < theDR_Veto_HO) {
0330         LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto HO tower: Calo deltaR= " << deltarHOcal;
0331         LogDebug("RecoMuon|CaloExtractorByAssociator")
0332             << " >>> Calo eta phi ethcal: " << calCPtr->eta() << " " << calCPtr->phi() << " " << ethcal;
0333         vetoTowerHOCal = true;
0334       }
0335 
0336       //! and now pitch those in the crossed list
0337       if (!(vetoTowerHOCal && vetoTowerHcal && vetoTowerEcal)) {
0338         for (unsigned int iH = 0; iH < mInfo.crossedTowerIds.size(); ++iH) {
0339           if (mInfo.crossedTowerIds[iH].rawId() == calCPtr->id().rawId()) {
0340             vetoTowerEcal = true;
0341             vetoTowerHcal = true;
0342             vetoTowerHOCal = true;
0343             break;
0344           }
0345         }
0346       }
0347 
0348       if (deltar0 > theDR_Max && !(vetoTowerEcal || vetoTowerHcal || vetoTowerHOCal))
0349         continue;
0350 
0351       reco::isodeposit::Direction towerDir(calCPtr->eta(), calCPtr->phi());
0352       //! add the Et of the tower to deposits if it's not a vetoed; put into muonEnergy otherwise
0353       if (doEcal) {
0354         if (vetoTowerEcal)
0355           depEcal.addCandEnergy(etecal);
0356         else if (deltar0 <= theDR_Max)
0357           depEcal.addDeposit(towerDir, etecal);
0358       }
0359       if (doHcal) {
0360         if (vetoTowerHcal)
0361           depHcal.addCandEnergy(ethcal);
0362         else if (deltar0 <= theDR_Max)
0363           depHcal.addDeposit(towerDir, ethcal);
0364       }
0365       if (doHOcal) {
0366         if (vetoTowerHOCal)
0367           depHOcal.addCandEnergy(ethocal);
0368         else if (deltar0 <= theDR_Max)
0369           depHOcal.addDeposit(towerDir, ethocal);
0370       }
0371     }
0372   }
0373 
0374   std::vector<IsoDeposit> resultDeps;
0375   resultDeps.push_back(depEcal);
0376   resultDeps.push_back(depHcal);
0377   resultDeps.push_back(depHOcal);
0378 
0379   return resultDeps;
0380 }
0381 
0382 double CaloExtractorByAssociator::noiseEcal(const CaloTower& tower) const {
0383   double noise = theNoiseTow_EB;
0384   double eta = tower.eta();
0385   if (fabs(eta) > 1.479)
0386     noise = theNoiseTow_EE;
0387   return noise;
0388 }
0389 
0390 double CaloExtractorByAssociator::noiseHcal(const CaloTower& tower) const {
0391   double noise = fabs(tower.eta()) > 1.479 ? theNoise_HE : theNoise_HB;
0392   return noise;
0393 }
0394 
0395 double CaloExtractorByAssociator::noiseHOcal(const CaloTower& tower) const {
0396   double noise = theNoise_HO;
0397   return noise;
0398 }
0399 
0400 double CaloExtractorByAssociator::noiseRecHit(const DetId& detId) const {
0401   double noise = 100;
0402   DetId::Detector det = detId.det();
0403   if (det == DetId::Ecal) {
0404     EcalSubdetector subDet = (EcalSubdetector)(detId.subdetId());
0405     if (subDet == EcalBarrel) {
0406       noise = theNoise_EB;
0407     } else if (subDet == EcalEndcap) {
0408       noise = theNoise_EE;
0409     }
0410   } else if (det == DetId::Hcal) {
0411     HcalSubdetector subDet = (HcalSubdetector)(detId.subdetId());
0412     if (subDet == HcalBarrel) {
0413       noise = theNoise_HB;
0414     } else if (subDet == HcalEndcap) {
0415       noise = theNoise_HE;
0416     } else if (subDet == HcalOuter) {
0417       noise = theNoise_HO;
0418     }
0419   }
0420   return noise;
0421 }
[ 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!