Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:04 

0001 #include "CaloExtractor.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/Math/interface/deltaR.h"
0008 #include "DataFormats/Math/interface/normalizedPhi.h"
0009 
0010 using namespace edm;
0011 using namespace std;
0012 using namespace reco;
0013 using namespace muonisolation;
0014 using reco::isodeposit::Direction;
0015 
0016 CaloExtractor::CaloExtractor(const ParameterSet& par, edm::ConsumesCollector&& iC)
0017     : theCaloTowerCollectionToken(
0018           iC.consumes<CaloTowerCollection>(par.getParameter<edm::InputTag>("CaloTowerCollectionLabel"))),
0019       theDepositLabel(par.getUntrackedParameter<string>("DepositLabel")),
0020       theCaloGeomToken(iC.esConsumes()),
0021       theFieldToken(iC.esConsumes()),
0022       theWeight_E(par.getParameter<double>("Weight_E")),
0023       theWeight_H(par.getParameter<double>("Weight_H")),
0024       theThreshold_E(par.getParameter<double>("Threshold_E")),
0025       theThreshold_H(par.getParameter<double>("Threshold_H")),
0026       theDR_Veto_E(par.getParameter<double>("DR_Veto_E")),
0027       theDR_Veto_H(par.getParameter<double>("DR_Veto_H")),
0028       theDR_Max(par.getParameter<double>("DR_Max")),
0029       vertexConstraintFlag_XY(par.getParameter<bool>("Vertex_Constraint_XY")),
0030       vertexConstraintFlag_Z(par.getParameter<bool>("Vertex_Constraint_Z")) {}
0031 
0032 void CaloExtractor::fillVetos(const edm::Event& event,
0033                               const edm::EventSetup& eventSetup,
0034                               const TrackCollection& muons) {
0035   theVetoCollection.clear();
0036 
0037   Handle<CaloTowerCollection> towers;
0038   event.getByToken(theCaloTowerCollectionToken, towers);
0039 
0040   auto const& caloGeom = eventSetup.getData(theCaloGeomToken);
0041 
0042   auto const& bField = eventSetup.getData(theFieldToken);
0043   double bz = bField.inInverseGeV(GlobalPoint(0., 0., 0.)).z();
0044 
0045   TrackCollection::const_iterator mu;
0046   TrackCollection::const_iterator muEnd(muons.end());
0047 
0048   CaloTowerCollection::const_iterator cal;
0049   CaloTowerCollection::const_iterator calEnd(towers->end());
0050 
0051   for (mu = muons.begin(); mu != muEnd; ++mu) {
0052     for (cal = towers->begin(); cal != calEnd; ++cal) {
0053       //! make this abit faster
0054       double dEta = fabs(mu->eta() - cal->eta());
0055       if (fabs(dEta) > theDR_Max)
0056         continue;
0057 
0058       double deltar0 = reco::deltaR(*mu, *cal);
0059       if (deltar0 > theDR_Max)
0060         continue;
0061 
0062       double etecal = cal->emEt();
0063       double eecal = cal->emEnergy();
0064       bool doEcal = theWeight_E > 0 && etecal > theThreshold_E && eecal > 3 * noiseEcal(*cal);
0065       double ethcal = cal->hadEt();
0066       double ehcal = cal->hadEnergy();
0067       bool doHcal = theWeight_H > 0 && ethcal > theThreshold_H && ehcal > 3 * noiseHcal(*cal);
0068       if ((!doEcal) && (!doHcal))
0069         continue;
0070 
0071       DetId calId = cal->id();
0072       GlobalPoint endpos = caloGeom.getPosition(calId);
0073       GlobalPoint muatcal = MuonAtCaloPosition(*mu, bz, endpos, vertexConstraintFlag_XY, vertexConstraintFlag_Z);
0074       double deltar = reco::deltaR(muatcal, endpos);
0075 
0076       if (doEcal) {
0077         if (deltar < theDR_Veto_E)
0078           theVetoCollection.push_back(calId);
0079       } else {
0080         if (deltar < theDR_Veto_H)
0081           theVetoCollection.push_back(calId);
0082       }
0083     }
0084   }
0085 }
0086 
0087 IsoDeposit CaloExtractor::deposit(const Event& event, const EventSetup& eventSetup, const Track& muon) const {
0088   IsoDeposit dep(muon.eta(), muon.phi());
0089   LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
0090       << " >>> Muon: pt " << muon.pt() << " eta " << muon.eta() << " phi " << muon.phi();
0091 
0092   Handle<CaloTowerCollection> towers;
0093   event.getByToken(theCaloTowerCollectionToken, towers);
0094 
0095   auto const& caloGeom = eventSetup.getData(theCaloGeomToken);
0096 
0097   auto const& bField = eventSetup.getData(theFieldToken);
0098   double bz = bField.inInverseGeV(GlobalPoint(0., 0., 0.)).z();
0099 
0100   CaloTowerCollection::const_iterator cal;
0101   CaloTowerCollection::const_iterator calEnd(towers->end());
0102   for (cal = towers->begin(); cal != calEnd; ++cal) {
0103     //! make this abit faster
0104     double dEta = fabs(muon.eta() - cal->eta());
0105     if (fabs(dEta) > theDR_Max)
0106       continue;
0107 
0108     double deltar0 = reco::deltaR(muon, *cal);
0109     if (deltar0 > theDR_Max)
0110       continue;
0111 
0112     double etecal = cal->emEt();
0113     double eecal = cal->emEnergy();
0114     bool doEcal = theWeight_E > 0 && etecal > theThreshold_E && eecal > 3 * noiseEcal(*cal);
0115     double ethcal = cal->hadEt();
0116     double ehcal = cal->hadEnergy();
0117     bool doHcal = theWeight_H > 0 && ethcal > theThreshold_H && ehcal > 3 * noiseHcal(*cal);
0118     if ((!doEcal) && (!doHcal))
0119       continue;
0120 
0121     DetId calId = cal->id();
0122     GlobalPoint endpos = caloGeom.getPosition(calId);
0123     GlobalPoint muatcal = MuonAtCaloPosition(muon, bz, endpos, vertexConstraintFlag_XY, vertexConstraintFlag_Z);
0124     double deltar = reco::deltaR(muatcal, endpos);
0125 
0126     if (deltar < theDR_Veto_H) {
0127       dep.setVeto(IsoDeposit::Veto(reco::isodeposit::Direction(muatcal.eta(), muatcal.phi()), theDR_Veto_H));
0128     }
0129 
0130     if (doEcal) {
0131       if (deltar < theDR_Veto_E) {
0132         double calodep = theWeight_E * etecal;
0133         if (doHcal)
0134           calodep += theWeight_H * ethcal;
0135         dep.addCandEnergy(calodep);
0136         LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
0137             << " >>> Calo deposit inside veto (with ECAL): deltar " << deltar << " calodep " << calodep << " ecaldep "
0138             << etecal << " hcaldep " << ethcal << " eta " << cal->eta() << " phi " << cal->phi();
0139         continue;
0140       }
0141     } else {
0142       if (deltar < theDR_Veto_H) {
0143         dep.addCandEnergy(theWeight_H * ethcal);
0144         LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
0145             << " >>> Calo deposit inside veto (no ECAL): deltar " << deltar << " calodep " << theWeight_H * ethcal
0146             << " eta " << cal->eta() << " phi " << cal->phi();
0147         continue;
0148       }
0149     }
0150 
0151     if (std::find(theVetoCollection.begin(), theVetoCollection.end(), calId) != theVetoCollection.end()) {
0152       LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
0153           << " >>> Deposits belongs to other track: deltar, etecal, ethcal= " << deltar << ", " << etecal << ", "
0154           << ethcal;
0155       continue;
0156     }
0157 
0158     if (doEcal) {
0159       if (deltar > theDR_Veto_E) {
0160         double calodep = theWeight_E * etecal;
0161         if (doHcal)
0162           calodep += theWeight_H * ethcal;
0163         dep.addDeposit(reco::isodeposit::Direction(endpos.eta(), endpos.phi()), calodep);
0164         LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
0165             << " >>> Calo deposit (with ECAL): deltar " << deltar << " calodep " << calodep << " ecaldep " << etecal
0166             << " hcaldep " << ethcal << " eta " << cal->eta() << " phi " << cal->phi();
0167       }
0168     } else {
0169       if (deltar > theDR_Veto_H) {
0170         dep.addDeposit(reco::isodeposit::Direction(endpos.eta(), endpos.phi()), theWeight_H * ethcal);
0171         LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
0172             << " >>> Calo deposit (no ECAL): deltar " << deltar << " calodep " << theWeight_H * ethcal << " eta "
0173             << cal->eta() << " phi " << cal->phi();
0174       }
0175     }
0176   }
0177 
0178   return dep;
0179 }
0180 
0181 GlobalPoint CaloExtractor::MuonAtCaloPosition(
0182     const Track& muon, const double bz, const GlobalPoint& endpos, bool fixVxy, bool fixVz) {
0183   double qoverp = muon.qoverp();
0184   double cur = bz * muon.charge() / muon.pt();
0185   double phi0 = muon.phi();
0186   double dca = muon.dxy();
0187   double theta = muon.theta();
0188   double dz = muon.dz();
0189 
0190   //LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
0191   //<< " Pt(GeV): " <<  muon.pt()
0192   //<< ", phi0 " <<  muon.phi0()
0193   //<< ", eta " <<  muon.eta();
0194   //LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
0195   //<< " d0 " <<  muon.d0()
0196   //<< ", dz " <<  muon.dz();
0197   //LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
0198   //<< " rhocal " <<  endpos.perp()
0199   //<< ", zcal " <<  endpos.z();
0200 
0201   if (fixVxy && fixVz) {
0202     // Note that here we assume no correlation between XY and Z projections
0203     // This should be a reasonable approximation for our purposes
0204     double errd02 = muon.covariance(muon.i_dxy, muon.i_dxy);
0205     if (pow(muon.dxy(), 2) < 4 * errd02) {
0206       phi0 -= muon.dxy() * muon.covariance(muon.i_dxy, muon.i_phi) / errd02;
0207       cur -= muon.dxy() * muon.covariance(muon.i_dxy, muon.i_qoverp) / errd02 * (cur / qoverp);
0208       dca = 0;
0209     }
0210     double errdsz2 = muon.covariance(muon.i_dsz, muon.i_dsz);
0211     if (pow(muon.dsz(), 2) < 4 * errdsz2) {
0212       theta += muon.dsz() * muon.covariance(muon.i_dsz, muon.i_lambda) / errdsz2;
0213       dz = 0;
0214     }
0215   } else if (fixVxy) {
0216     double errd02 = muon.covariance(muon.i_dxy, muon.i_dxy);
0217     if (pow(muon.dxy(), 2) < 4 * errd02) {
0218       phi0 -= muon.dxy() * muon.covariance(muon.i_dxy, muon.i_phi) / errd02;
0219       cur -= muon.dxy() * muon.covariance(muon.i_dxy, muon.i_qoverp) / errd02 * (cur / qoverp);
0220       theta += muon.dxy() * muon.covariance(muon.i_dxy, muon.i_lambda) / errd02;
0221       dz -= muon.dxy() * muon.covariance(muon.i_dxy, muon.i_dsz) / errd02 * muon.p() / muon.pt();
0222       dca = 0;
0223     }
0224   } else if (fixVz) {
0225     double errdsz2 = muon.covariance(muon.i_dsz, muon.i_dsz);
0226     if (pow(muon.dsz(), 2) < 4 * errdsz2) {
0227       theta += muon.dsz() * muon.covariance(muon.i_dsz, muon.i_lambda) / errdsz2;
0228       phi0 -= muon.dsz() * muon.covariance(muon.i_dsz, muon.i_phi) / errdsz2;
0229       cur -= muon.dsz() * muon.covariance(muon.i_dsz, muon.i_qoverp) / errdsz2 * (cur / qoverp);
0230       dca -= muon.dsz() * muon.covariance(muon.i_dsz, muon.i_dxy) / errdsz2;
0231       dz = 0;
0232     }
0233   }
0234 
0235   double sphi0 = sin(phi0);
0236   double cphi0 = cos(phi0);
0237 
0238   double xsin = endpos.x() * sphi0 - endpos.y() * cphi0;
0239   double xcos = endpos.x() * cphi0 + endpos.y() * sphi0;
0240   double fcdca = fabs(1 - cur * dca);
0241   double phif = atan2(fcdca * sphi0 - cur * endpos.x(), fcdca * cphi0 + cur * endpos.y());
0242   double tphif2 = tan(0.5 * (phif - phi0));
0243   double dcaf = dca + xsin + xcos * tphif2;
0244 
0245   double x = endpos.x() - dcaf * sin(phif);
0246   double y = endpos.y() + dcaf * cos(phif);
0247 
0248   double deltas = (x - muon.vx()) * cphi0 + (y - muon.vy()) * sphi0;
0249   double deltaphi = normalizedPhi(phif - phi0);
0250   if (deltaphi != 0)
0251     deltas = deltas * deltaphi / sin(deltaphi);
0252 
0253   double z = dz;
0254   double tantheta = tan(theta);
0255   if (tantheta != 0) {
0256     z += deltas / tan(theta);
0257   } else {
0258     z = endpos.z();
0259   }
0260 
0261   return GlobalPoint(x, y, z);
0262 }
0263 
0264 double CaloExtractor::noiseEcal(const CaloTower& tower) const {
0265   double noise = 0.04;
0266   double eta = tower.eta();
0267   if (fabs(eta) > 1.479)
0268     noise = 0.15;
0269   return noise;
0270 }
0271 
0272 double CaloExtractor::noiseHcal(const CaloTower& tower) const {
0273   double noise = 0.2;
0274   return noise;
0275 }
[ 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!