CMS 3D CMS Logo
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
CaloExtractor.cc
Go to the documentation of this file.
1 #include "CaloExtractor.h"
2 
3 #include "DataFormats/Common/interface/Handle.h"
4 #include "FWCore/Framework/interface/ESHandle.h"
5 #include "FWCore/MessageLogger/interface/MessageLogger.h"
6 #include "FWCore/ParameterSet/interface/ParameterSet.h"
7 
8 #include "DataFormats/Math/interface/deltaR.h"
9 #include "DataFormats/Math/interface/normalizedPhi.h"
10 
11 using namespace edm;
12 using namespace std;
13 using namespace reco;
14 using namespace muonisolation;
15 using reco::isodeposit::Direction;
16 
17 CaloExtractor::CaloExtractor(const ParameterSet& par, edm::ConsumesCollector&& iC)
18  : theCaloTowerCollectionToken(
19  iC.consumes<CaloTowerCollection>(par.getParameter<edm::InputTag>("CaloTowerCollectionLabel"))),
20  theDepositLabel(par.getUntrackedParameter<string>("DepositLabel")),
21  theCaloGeomToken(iC.esConsumes()),
22  theFieldToken(iC.esConsumes()),
23  theWeight_E(par.getParameter<double>("Weight_E")),
24  theWeight_H(par.getParameter<double>("Weight_H")),
25  theThreshold_E(par.getParameter<double>("Threshold_E")),
26  theThreshold_H(par.getParameter<double>("Threshold_H")),
27  theDR_Veto_E(par.getParameter<double>("DR_Veto_E")),
28  theDR_Veto_H(par.getParameter<double>("DR_Veto_H")),
29  theDR_Max(par.getParameter<double>("DR_Max")),
30  vertexConstraintFlag_XY(par.getParameter<bool>("Vertex_Constraint_XY")),
31  vertexConstraintFlag_Z(par.getParameter<bool>("Vertex_Constraint_Z")) {}
32 
33 void CaloExtractor::fillVetos(const edm::Event& event,
34  const edm::EventSetup& eventSetup,
35  const TrackCollection& muons) {
36  theVetoCollection.clear();
37 
38  Handle<CaloTowerCollection> towers;
39  event.getByToken(theCaloTowerCollectionToken, towers);
40 
41  auto const& caloGeom = eventSetup.getData(theCaloGeomToken);
42 
43  auto const& bField = eventSetup.getData(theFieldToken);
44  double bz = bField.inInverseGeV(GlobalPoint(0., 0., 0.)).z();
45 
46  TrackCollection::const_iterator mu;
47  TrackCollection::const_iterator muEnd(muons.end());
48 
49  CaloTowerCollection::const_iterator cal;
50  CaloTowerCollection::const_iterator calEnd(towers->end());
51 
52  for (mu = muons.begin(); mu != muEnd; ++mu) {
53  for (cal = towers->begin(); cal != calEnd; ++cal) {
55  double dEta = fabs(mu->eta() - cal->eta());
56  if (fabs(dEta) > theDR_Max)
57  continue;
58 
59  double deltar0 = reco::deltaR(*mu, *cal);
60  if (deltar0 > theDR_Max)
61  continue;
62 
63  double etecal = cal->emEt();
64  double eecal = cal->emEnergy();
65  bool doEcal = theWeight_E > 0 && etecal > theThreshold_E && eecal > 3 * noiseEcal(*cal);
66  double ethcal = cal->hadEt();
67  double ehcal = cal->hadEnergy();
68  bool doHcal = theWeight_H > 0 && ethcal > theThreshold_H && ehcal > 3 * noiseHcal(*cal);
69  if ((!doEcal) && (!doHcal))
70  continue;
71 
72  DetId calId = cal->id();
73  GlobalPoint endpos = caloGeom.getPosition(calId);
74  GlobalPoint muatcal = MuonAtCaloPosition(*mu, bz, endpos, vertexConstraintFlag_XY, vertexConstraintFlag_Z);
75  double deltar = reco::deltaR(muatcal, endpos);
76 
77  if (doEcal) {
78  if (deltar < theDR_Veto_E)
79  theVetoCollection.push_back(calId);
80  } else {
81  if (deltar < theDR_Veto_H)
82  theVetoCollection.push_back(calId);
83  }
84  }
85  }
86 }
87 
88 IsoDeposit CaloExtractor::deposit(const Event& event, const EventSetup& eventSetup, const Track& muon) const {
89  IsoDeposit dep(muon.eta(), muon.phi());
90  LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
91  << " >>> Muon: pt " << muon.pt() << " eta " << muon.eta() << " phi " << muon.phi();
92 
93  Handle<CaloTowerCollection> towers;
94  event.getByToken(theCaloTowerCollectionToken, towers);
95 
96  auto const& caloGeom = eventSetup.getData(theCaloGeomToken);
97 
98  auto const& bField = eventSetup.getData(theFieldToken);
99  double bz = bField.inInverseGeV(GlobalPoint(0., 0., 0.)).z();
100 
101  CaloTowerCollection::const_iterator cal;
102  CaloTowerCollection::const_iterator calEnd(towers->end());
103  for (cal = towers->begin(); cal != calEnd; ++cal) {
105  double dEta = fabs(muon.eta() - cal->eta());
106  if (fabs(dEta) > theDR_Max)
107  continue;
108 
109  double deltar0 = reco::deltaR(muon, *cal);
110  if (deltar0 > theDR_Max)
111  continue;
112 
113  double etecal = cal->emEt();
114  double eecal = cal->emEnergy();
115  bool doEcal = theWeight_E > 0 && etecal > theThreshold_E && eecal > 3 * noiseEcal(*cal);
116  double ethcal = cal->hadEt();
117  double ehcal = cal->hadEnergy();
118  bool doHcal = theWeight_H > 0 && ethcal > theThreshold_H && ehcal > 3 * noiseHcal(*cal);
119  if ((!doEcal) && (!doHcal))
120  continue;
121 
122  DetId calId = cal->id();
123  GlobalPoint endpos = caloGeom.getPosition(calId);
124  GlobalPoint muatcal = MuonAtCaloPosition(muon, bz, endpos, vertexConstraintFlag_XY, vertexConstraintFlag_Z);
125  double deltar = reco::deltaR(muatcal, endpos);
126 
127  if (deltar < theDR_Veto_H) {
128  dep.setVeto(IsoDeposit::Veto(reco::isodeposit::Direction(muatcal.eta(), muatcal.phi()), theDR_Veto_H));
129  }
130 
131  if (doEcal) {
132  if (deltar < theDR_Veto_E) {
133  double calodep = theWeight_E * etecal;
134  if (doHcal)
135  calodep += theWeight_H * ethcal;
136  dep.addCandEnergy(calodep);
137  LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
138  << " >>> Calo deposit inside veto (with ECAL): deltar " << deltar << " calodep " << calodep << " ecaldep "
139  << etecal << " hcaldep " << ethcal << " eta " << cal->eta() << " phi " << cal->phi();
140  continue;
141  }
142  } else {
143  if (deltar < theDR_Veto_H) {
144  dep.addCandEnergy(theWeight_H * ethcal);
145  LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
146  << " >>> Calo deposit inside veto (no ECAL): deltar " << deltar << " calodep " << theWeight_H * ethcal
147  << " eta " << cal->eta() << " phi " << cal->phi();
148  continue;
149  }
150  }
151 
152  if (std::find(theVetoCollection.begin(), theVetoCollection.end(), calId) != theVetoCollection.end()) {
153  LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
154  << " >>> Deposits belongs to other track: deltar, etecal, ethcal= " << deltar << ", " << etecal << ", "
155  << ethcal;
156  continue;
157  }
158 
159  if (doEcal) {
160  if (deltar > theDR_Veto_E) {
161  double calodep = theWeight_E * etecal;
162  if (doHcal)
163  calodep += theWeight_H * ethcal;
164  dep.addDeposit(reco::isodeposit::Direction(endpos.eta(), endpos.phi()), calodep);
165  LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
166  << " >>> Calo deposit (with ECAL): deltar " << deltar << " calodep " << calodep << " ecaldep " << etecal
167  << " hcaldep " << ethcal << " eta " << cal->eta() << " phi " << cal->phi();
168  }
169  } else {
170  if (deltar > theDR_Veto_H) {
171  dep.addDeposit(reco::isodeposit::Direction(endpos.eta(), endpos.phi()), theWeight_H * ethcal);
172  LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
173  << " >>> Calo deposit (no ECAL): deltar " << deltar << " calodep " << theWeight_H * ethcal << " eta "
174  << cal->eta() << " phi " << cal->phi();
175  }
176  }
177  }
178 
179  return dep;
180 }
181 
182 GlobalPoint CaloExtractor::MuonAtCaloPosition(
183  const Track& muon, const double bz, const GlobalPoint& endpos, bool fixVxy, bool fixVz) {
184  double qoverp = muon.qoverp();
185  double cur = bz * muon.charge() / muon.pt();
186  double phi0 = muon.phi();
187  double dca = muon.dxy();
188  double theta = muon.theta();
189  double dz = muon.dz();
190 
191  //LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
192  //<< " Pt(GeV): " << muon.pt()
193  //<< ", phi0 " << muon.phi0()
194  //<< ", eta " << muon.eta();
195  //LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
196  //<< " d0 " << muon.d0()
197  //<< ", dz " << muon.dz();
198  //LogDebug("Muon|RecoMuon|L2MuonIsolationProducer")
199  //<< " rhocal " << endpos.perp()
200  //<< ", zcal " << endpos.z();
201 
202  if (fixVxy && fixVz) {
203  // Note that here we assume no correlation between XY and Z projections
204  // This should be a reasonable approximation for our purposes
205  double errd02 = muon.covariance(muon.i_dxy, muon.i_dxy);
206  if (pow(muon.dxy(), 2) < 4 * errd02) {
207  phi0 -= muon.dxy() * muon.covariance(muon.i_dxy, muon.i_phi) / errd02;
208  cur -= muon.dxy() * muon.covariance(muon.i_dxy, muon.i_qoverp) / errd02 * (cur / qoverp);
209  dca = 0;
210  }
211  double errdsz2 = muon.covariance(muon.i_dsz, muon.i_dsz);
212  if (pow(muon.dsz(), 2) < 4 * errdsz2) {
213  theta += muon.dsz() * muon.covariance(muon.i_dsz, muon.i_lambda) / errdsz2;
214  dz = 0;
215  }
216  } else if (fixVxy) {
217  double errd02 = muon.covariance(muon.i_dxy, muon.i_dxy);
218  if (pow(muon.dxy(), 2) < 4 * errd02) {
219  phi0 -= muon.dxy() * muon.covariance(muon.i_dxy, muon.i_phi) / errd02;
220  cur -= muon.dxy() * muon.covariance(muon.i_dxy, muon.i_qoverp) / errd02 * (cur / qoverp);
221  theta += muon.dxy() * muon.covariance(muon.i_dxy, muon.i_lambda) / errd02;
222  dz -= muon.dxy() * muon.covariance(muon.i_dxy, muon.i_dsz) / errd02 * muon.p() / muon.pt();
223  dca = 0;
224  }
225  } else if (fixVz) {
226  double errdsz2 = muon.covariance(muon.i_dsz, muon.i_dsz);
227  if (pow(muon.dsz(), 2) < 4 * errdsz2) {
228  theta += muon.dsz() * muon.covariance(muon.i_dsz, muon.i_lambda) / errdsz2;
229  phi0 -= muon.dsz() * muon.covariance(muon.i_dsz, muon.i_phi) / errdsz2;
230  cur -= muon.dsz() * muon.covariance(muon.i_dsz, muon.i_qoverp) / errdsz2 * (cur / qoverp);
231  dca -= muon.dsz() * muon.covariance(muon.i_dsz, muon.i_dxy) / errdsz2;
232  dz = 0;
233  }
234  }
235 
236  double sphi0 = sin(phi0);
237  double cphi0 = cos(phi0);
238 
239  double xsin = endpos.x() * sphi0 - endpos.y() * cphi0;
240  double xcos = endpos.x() * cphi0 + endpos.y() * sphi0;
241  double fcdca = fabs(1 - cur * dca);
242  double phif = atan2(fcdca * sphi0 - cur * endpos.x(), fcdca * cphi0 + cur * endpos.y());
243  double tphif2 = tan(0.5 * (phif - phi0));
244  double dcaf = dca + xsin + xcos * tphif2;
245 
246  double x = endpos.x() - dcaf * sin(phif);
247  double y = endpos.y() + dcaf * cos(phif);
248 
249  double deltas = (x - muon.vx()) * cphi0 + (y - muon.vy()) * sphi0;
250  double deltaphi = normalizedPhi(phif - phi0);
251  if (deltaphi != 0)
252  deltas = deltas * deltaphi / sin(deltaphi);
253 
254  double z = dz;
255  double tantheta = tan(theta);
256  if (tantheta != 0) {
257  z += deltas / tan(theta);
258  } else {
259  z = endpos.z();
260  }
261 
262  return GlobalPoint(x, y, z);
263 }
264 
265 double CaloExtractor::noiseEcal(const CaloTower& tower) const {
266  double noise = 0.04;
267  double eta = tower.eta();
268  if (fabs(eta) > 1.479)
269  noise = 0.15;
270  return noise;
271 }
272 
273 double CaloExtractor::noiseHcal(const CaloTower& tower) const {
274  double noise = 0.2;
275  return noise;
276 }
PDWG_BPHSkim_cff.muons
muons
Definition: PDWG_BPHSkim_cff.py:47
HLT_FULL_cff.towers
towers
Definition: HLT_FULL_cff.py:36379
muonisolation::CaloExtractor::noiseEcal
double noiseEcal(const CaloTower &tower) const
Definition: CaloExtractor.cc:265
Handle.h
electrons_cff.bool
bool
Definition: electrons_cff.py:366
edm::SortedCollection< CaloTower >::const_iterator
std::vector< CaloTower >::const_iterator const_iterator
Definition: SortedCollection.h:80
MessageLogger.h
ESHandle.h
muon
Definition: MuonCocktails.h:17
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
amptDefaultParameters_cff.mu
mu
Definition: amptDefaultParameters_cff.py:16
edm
HLT enums.
Definition: AlignableModifier.h:19
muonisolation::CaloExtractor::theCaloGeomToken
const edm::ESGetToken< CaloGeometry, CaloGeometryRecord > theCaloGeomToken
Definition: CaloExtractor.h:48
CaloExtractor.h
muonisolation::CaloExtractor::theThreshold_H
double theThreshold_H
Definition: CaloExtractor.h:55
hgcalTowerProducer_cfi.tower
tower
Definition: hgcalTowerProducer_cfi.py:4
edm::SortedCollection< CaloTower >
deltar
Definition: deltar.py:1
reco::IsoDeposit::Veto
Definition: IsoDeposit.h:59
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:46
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
edm::Handle
Definition: AssociativeIterator.h:50
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
DetId
Definition: DetId.h:17
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
muonisolation::CaloExtractor::theThreshold_E
double theThreshold_E
Definition: CaloExtractor.h:54
muonisolation::CaloExtractor::theDR_Max
double theDR_Max
Definition: CaloExtractor.h:58
muonisolation::CaloExtractor::theFieldToken
const edm::ESGetToken< MagneticField, IdealMagneticFieldRecord > theFieldToken
Definition: CaloExtractor.h:49
PVValHelper::eta
Definition: PVValidationHelpers.h:70
muonisolation::CaloExtractor::theVetoCollection
std::vector< DetId > theVetoCollection
Definition: CaloExtractor.h:63
muonisolation::CaloExtractor::theDR_Veto_E
double theDR_Veto_E
Definition: CaloExtractor.h:56
reco::Track
Definition: Track.h:27
normalizedPhi
constexpr T normalizedPhi(T phi)
Definition: normalizedPhi.h:8
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
GlobalPoint
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
Point3DBase< float, GlobalTag >
muonisolation::CaloExtractor::theWeight_E
double theWeight_E
Definition: CaloExtractor.h:52
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:233
edm::ParameterSet
Definition: ParameterSet.h:47
deltaR.h
muonisolation::CaloExtractor::theWeight_H
double theWeight_H
Definition: CaloExtractor.h:53
muonisolation::CaloExtractor::deposit
reco::IsoDeposit deposit(const edm::Event &ev, const edm::EventSetup &evSetup, const reco::Track &track) const override
Definition: CaloExtractor.cc:88
muonisolation::CaloExtractor::noiseHcal
double noiseHcal(const CaloTower &tower) const
Definition: CaloExtractor.cc:273
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:73
muonisolation::CaloExtractor::theCaloTowerCollectionToken
edm::EDGetTokenT< CaloTowerCollection > theCaloTowerCollectionToken
Definition: CaloExtractor.h:43
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
normalizedPhi.h
CaloTower
Definition: CaloTower.h:26
edm::EventSetup
Definition: EventSetup.h:58
gpuVertexFinder::noise
__shared__ int noise
Definition: gpuFitVertices.h:50
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
muonisolation::CaloExtractor::vertexConstraintFlag_XY
bool vertexConstraintFlag_XY
Definition: CaloExtractor.h:59
edm::EventSetup::getData
bool getData(T &iHolder) const
Definition: EventSetup.h:127
muonisolation
Definition: CandViewExtractor.h:16
std
Definition: JetResolutionObject.h:76
muonisolation::CaloExtractor::vertexConstraintFlag_Z
bool vertexConstraintFlag_Z
Definition: CaloExtractor.h:60
Calorimetry_cff.bField
bField
Definition: Calorimetry_cff.py:284
muonisolation::CaloExtractor::theDR_Veto_H
double theDR_Veto_H
Definition: CaloExtractor.h:57
PVValHelper::dz
Definition: PVValidationHelpers.h:51
HLT_FULL_cff.dEta
dEta
Definition: HLT_FULL_cff.py:13694
reco::isodeposit::Direction
Definition: IsoDepositDirection.h:19
muonisolation::CaloExtractor::fillVetos
void fillVetos(const edm::Event &ev, const edm::EventSetup &evSetup, const reco::TrackCollection &tracks) override
Definition: CaloExtractor.cc:33
reco::IsoDeposit
Definition: IsoDeposit.h:49
reco::deltaR
constexpr auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:30
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
ParameterSet.h
event
Definition: event.py:1
DeDxTools::esConsumes
ESGetTokenH3DDVariant esConsumes(std::string const &Reccord, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
edm::Event
Definition: Event.h:73
muonisolation::CaloExtractor::MuonAtCaloPosition
static GlobalPoint MuonAtCaloPosition(const reco::Track &muon, const double bz, const GlobalPoint &endpos, bool fixVxy=false, bool fixVz=false)
Extrapolate muons to calorimeter-object positions.
Definition: CaloExtractor.cc:182
edm::InputTag
Definition: InputTag.h:15
reco::TrackCollection
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
edm::ConsumesCollector
Definition: ConsumesCollector.h:45
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
Generated for CMSSW Reference Manual by doxygen 1.8.16