CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CaloExtractorByAssociator.cc
Go to the documentation of this file.
2 
10 
15 
20 
22 
26 
28 //#include "CommonTools/Utils/interface/deltaR.h"
29 //#include "PhysicsTools/Utilities/interface/deltaR.h"
30 
31 using namespace edm;
32 using namespace std;
33 using namespace reco;
34 using namespace muonisolation;
36 
37 CaloExtractorByAssociator::CaloExtractorByAssociator(const ParameterSet& par) :
38  theUseRecHitsFlag(par.getParameter<bool>("UseRecHitsFlag")),
39  theDepositLabel(par.getUntrackedParameter<string>("DepositLabel")),
40  theDepositInstanceLabels(par.getParameter<std::vector<std::string> >("DepositInstanceLabels")),
41  thePropagatorName(par.getParameter<std::string>("PropagatorName")),
42  theThreshold_E(par.getParameter<double>("Threshold_E")),
43  theThreshold_H(par.getParameter<double>("Threshold_H")),
44  theThreshold_HO(par.getParameter<double>("Threshold_HO")),
45  theDR_Veto_E(par.getParameter<double>("DR_Veto_E")),
46  theDR_Veto_H(par.getParameter<double>("DR_Veto_H")),
47  theDR_Veto_HO(par.getParameter<double>("DR_Veto_HO")),
48  theCenterConeOnCalIntersection(par.getParameter<bool>("CenterConeOnCalIntersection")),
49  theDR_Max(par.getParameter<double>("DR_Max")),
50  theNoise_EB(par.getParameter<double>("Noise_EB")),
51  theNoise_EE(par.getParameter<double>("Noise_EE")),
52  theNoise_HB(par.getParameter<double>("Noise_HB")),
53  theNoise_HE(par.getParameter<double>("Noise_HE")),
54  theNoise_HO(par.getParameter<double>("Noise_HO")),
55  theNoiseTow_EB(par.getParameter<double>("NoiseTow_EB")),
56  theNoiseTow_EE(par.getParameter<double>("NoiseTow_EE")),
57  theService(0),
58  theAssociator(0),
59  thePrintTimeReport(par.getUntrackedParameter<bool>("PrintTimeReport"))
60 {
61  ParameterSet serviceParameters = par.getParameter<ParameterSet>("ServiceParameters");
62  theService = new MuonServiceProxy(serviceParameters);
63 
66 }
67 
71  if (theService) delete theService;
72  if (theAssociator) delete theAssociator;
73 }
74 
76 {
77 // LogWarning("CaloExtractorByAssociator")
78 // <<"fillVetos does nothing now: IsoDeposit provides enough functionality\n"
79 // <<"to remove a deposit at/around given (eta, phi)";
80 
81 }
82 
83 IsoDeposit CaloExtractorByAssociator::deposit( const Event & event, const EventSetup& eventSetup, const Track & muon) const
84 {
85  IsoDeposit::Direction muonDir(muon.eta(), muon.phi());
86  IsoDeposit dep(muonDir );
87 
88 // LogWarning("CaloExtractorByAssociator")
89 // <<"single deposit is not an option here\n"
90 // <<"use ::deposits --> extract all and reweight as necessary";
91 
92  return dep;
93 
94 }
95 
96 
98 std::vector<IsoDeposit> CaloExtractorByAssociator::deposits( const Event & event, const EventSetup& eventSetup, const Track & muon) const
99 {
100  theService->update(eventSetup);
102 
105  if (theDepositInstanceLabels.size() != 3){
106  LogError("MuonIsolation")<<"Configuration is inconsistent: Need 3 deposit instance labels";
107  }
108  if (! theDepositInstanceLabels[0].compare(0,1, std::string("e")) == 0
109  || ! theDepositInstanceLabels[1].compare(0,1, std::string("h")) == 0
110  || ! theDepositInstanceLabels[2].compare(0,2, std::string("ho")) == 0){
111  LogWarning("MuonIsolation")<<"Deposit instance labels do not look like (e*, h*, ho*):"
112  <<"proceed at your own risk. The extractor interprets lab0=from ecal; lab1=from hcal; lab2=from ho";
113  }
114 
115  typedef IsoDeposit::Veto Veto;
118  IsoDeposit::Direction muonDir(muon.eta(), muon.phi());
119 
120  IsoDeposit depEcal(muonDir);
121  IsoDeposit depHcal(muonDir);
122  IsoDeposit depHOcal(muonDir);
123 
125  eventSetup.get<IdealMagneticFieldRecord>().get(bField);
126 
127 
128  reco::TransientTrack tMuon(muon, &*bField);
129  FreeTrajectoryState iFTS = tMuon.initialFreeState();
130  TrackDetMatchInfo mInfo = theAssociator->associate(event, eventSetup, iFTS, *theAssociatorParameters);
131 
133  depEcal.setVeto(Veto(reco::isodeposit::Direction(mInfo.trkGlobPosAtEcal.eta(), mInfo.trkGlobPosAtEcal.phi()),
134  theDR_Veto_E));
135  depHcal.setVeto(Veto(reco::isodeposit::Direction(mInfo.trkGlobPosAtHcal.eta(), mInfo.trkGlobPosAtHcal.phi()),
136  theDR_Veto_H));
137  depHOcal.setVeto(Veto(reco::isodeposit::Direction(mInfo.trkGlobPosAtHO.eta(), mInfo.trkGlobPosAtHO.phi()),
138  theDR_Veto_HO));
139 
141  reco::isodeposit::Direction dirTmp = depEcal.veto().vetoDir;
142  double dRtmp = depEcal.veto().dR;
143  depEcal = IsoDeposit(dirTmp); depEcal.setVeto(Veto(dirTmp, dRtmp));
144 
145  dirTmp = depHcal.veto().vetoDir;
146  dRtmp = depHcal.veto().dR;
147  depHcal = IsoDeposit(dirTmp); depHcal.setVeto(Veto(dirTmp, dRtmp));
148 
149  dirTmp = depHOcal.veto().vetoDir;
150  dRtmp = depHOcal.veto().dR;
151  depHOcal = IsoDeposit(dirTmp); depHOcal.setVeto(Veto(dirTmp, dRtmp));
152  }
153 
154  if (theUseRecHitsFlag){
158  eventSetup.get<CaloGeometryRecord>().get(caloGeom);
159 
160  //Ecal
161  std::vector<const EcalRecHit*>::const_iterator eHitCI = mInfo.ecalRecHits.begin();
162  for (; eHitCI != mInfo.ecalRecHits.end(); ++eHitCI){
163  const EcalRecHit* eHitCPtr = *eHitCI;
164  GlobalPoint eHitPos = caloGeom->getPosition(eHitCPtr->detid());
165  double deltar0 = reco::deltaR(muon, eHitPos);
166  double cosTheta = 1./cosh(eHitPos.eta());
167  double energy = eHitCPtr->energy();
168  double et = energy*cosTheta;
169  if (deltar0 > theDR_Max
170  || ! (et > theThreshold_E && energy > 3*noiseRecHit(eHitCPtr->detid()))) continue;
171 
172  bool vetoHit = false;
173  double deltar = reco::deltaR(mInfo.trkGlobPosAtEcal, eHitPos);
175  if (deltar < theDR_Veto_E ){
176  LogDebug("RecoMuon|CaloExtractorByAssociator")
177  << " >>> Veto ECAL hit: Calo deltaR= " << deltar;
178  LogDebug("RecoMuon|CaloExtractorByAssociator")
179  << " >>> Calo eta phi ethcal: " << eHitPos.eta() << " " << eHitPos.phi() << " " << et;
180  vetoHit = true;
181  }
183  if (! vetoHit){
184  for (unsigned int iH = 0; iH< mInfo.crossedEcalIds.size() && ! vetoHit; ++iH){
185  if (mInfo.crossedEcalIds[iH].rawId() == eHitCPtr->detid().rawId()) vetoHit = true;
186  }
187  }
188 
189  if (vetoHit ){
190  depEcal.addCandEnergy(et);
191  } else {
192  depEcal.addDeposit(reco::isodeposit::Direction(eHitPos.eta(), eHitPos.phi()), et);
193  }
194  }
195 
196  //Hcal
197  std::vector<const HBHERecHit*>::const_iterator hHitCI = mInfo.hcalRecHits.begin();
198  for (; hHitCI != mInfo.hcalRecHits.end(); ++hHitCI){
199  const HBHERecHit* hHitCPtr = *hHitCI;
200  GlobalPoint hHitPos = caloGeom->getPosition(hHitCPtr->detid());
201  double deltar0 = reco::deltaR(muon, hHitPos);
202  double cosTheta = 1./cosh(hHitPos.eta());
203  double energy = hHitCPtr->energy();
204  double et = energy*cosTheta;
205  if (deltar0 > theDR_Max
206  || ! (et > theThreshold_H && energy > 3*noiseRecHit(hHitCPtr->detid()))) continue;
207 
208  bool vetoHit = false;
209  double deltar = reco::deltaR(mInfo.trkGlobPosAtHcal, hHitPos);
211  if (deltar < theDR_Veto_H ){
212  LogDebug("RecoMuon|CaloExtractorByAssociator")
213  << " >>> Veto HBHE hit: Calo deltaR= " << deltar;
214  LogDebug("RecoMuon|CaloExtractorByAssociator")
215  << " >>> Calo eta phi ethcal: " << hHitPos.eta() << " " << hHitPos.phi() << " " << et;
216  vetoHit = true;
217  }
219  if (! vetoHit){
220  for (unsigned int iH = 0; iH< mInfo.crossedHcalIds.size() && ! vetoHit; ++iH){
221  if (mInfo.crossedHcalIds[iH].rawId() == hHitCPtr->detid().rawId()) vetoHit = true;
222  }
223  }
224 
225  if (vetoHit ){
226  depHcal.addCandEnergy(et);
227  } else {
228  depHcal.addDeposit(reco::isodeposit::Direction(hHitPos.eta(), hHitPos.phi()), et);
229  }
230  }
231 
232  //HOcal
233  std::vector<const HORecHit*>::const_iterator hoHitCI = mInfo.hoRecHits.begin();
234  for (; hoHitCI != mInfo.hoRecHits.end(); ++hoHitCI){
235  const HORecHit* hoHitCPtr = *hoHitCI;
236  GlobalPoint hoHitPos = caloGeom->getPosition(hoHitCPtr->detid());
237  double deltar0 = reco::deltaR(muon, hoHitPos);
238  double cosTheta = 1./cosh(hoHitPos.eta());
239  double energy = hoHitCPtr->energy();
240  double et = energy*cosTheta;
241  if (deltar0 > theDR_Max
242  || ! (et > theThreshold_HO && energy > 3*noiseRecHit(hoHitCPtr->detid()))) continue;
243 
244  bool vetoHit = false;
245  double deltar = reco::deltaR(mInfo.trkGlobPosAtHO, hoHitPos);
247  if (deltar < theDR_Veto_HO ){
248  LogDebug("RecoMuon|CaloExtractorByAssociator")
249  << " >>> Veto HO hit: Calo deltaR= " << deltar;
250  LogDebug("RecoMuon|CaloExtractorByAssociator")
251  << " >>> Calo eta phi ethcal: " << hoHitPos.eta() << " " << hoHitPos.phi() << " " << et;
252  vetoHit = true;
253  }
255  if (! vetoHit){
256  for (unsigned int iH = 0; iH< mInfo.crossedHOIds.size() && ! vetoHit; ++iH){
257  if (mInfo.crossedHOIds[iH].rawId() == hoHitCPtr->detid().rawId()) vetoHit = true;
258  }
259  }
260 
261  if (vetoHit ){
262  depHOcal.addCandEnergy(et);
263  } else {
264  depHOcal.addDeposit(reco::isodeposit::Direction(hoHitPos.eta(), hoHitPos.phi()), et);
265  }
266  }
267 
268 
269  } else {
271  std::vector<const CaloTower*>::const_iterator calCI = mInfo.towers.begin();
272  for (; calCI != mInfo.towers.end(); ++calCI){
273  const CaloTower* calCPtr = *calCI;
274  double deltar0 = reco::deltaR(muon,*calCPtr);
275  if (deltar0>theDR_Max) continue;
276 
277  //even more copy-pasting .. need to refactor
278  double etecal = calCPtr->emEt();
279  double eecal = calCPtr->emEnergy();
280  bool doEcal = etecal>theThreshold_E && eecal>3*noiseEcal(*calCPtr);
281  double ethcal = calCPtr->hadEt();
282  double ehcal = calCPtr->hadEnergy();
283  bool doHcal = ethcal>theThreshold_H && ehcal>3*noiseHcal(*calCPtr);
284  double ethocal = calCPtr->outerEt();
285  double ehocal = calCPtr->outerEnergy();
286  bool doHOcal = ethocal>theThreshold_HO && ehocal>3*noiseHOcal(*calCPtr);
287  if ((!doEcal) && (!doHcal) && (!doHcal)) continue;
288 
289  bool vetoTowerEcal = false;
290  double deltarEcal = reco::deltaR(mInfo.trkGlobPosAtEcal, *calCPtr);
292  if (deltarEcal < theDR_Veto_E ){
293  LogDebug("RecoMuon|CaloExtractorByAssociator")
294  << " >>> Veto ecal tower: Calo deltaR= " << deltarEcal;
295  LogDebug("RecoMuon|CaloExtractorByAssociator")
296  << " >>> Calo eta phi ethcal: " << calCPtr->eta() << " " << calCPtr->phi() << " " << ethcal;
297  vetoTowerEcal = true;
298  }
299  bool vetoTowerHcal = false;
300  double deltarHcal = reco::deltaR(mInfo.trkGlobPosAtHcal, *calCPtr);
302  if (deltarHcal < theDR_Veto_H ){
303  LogDebug("RecoMuon|CaloExtractorByAssociator")
304  << " >>> Veto hcal tower: Calo deltaR= " << deltarHcal;
305  LogDebug("RecoMuon|CaloExtractorByAssociator")
306  << " >>> Calo eta phi ethcal: " << calCPtr->eta() << " " << calCPtr->phi() << " " << ethcal;
307  vetoTowerHcal = true;
308  }
309  bool vetoTowerHOCal = false;
310  double deltarHOcal = reco::deltaR(mInfo.trkGlobPosAtHO, *calCPtr);
312  if (deltarHOcal < theDR_Veto_HO ){
313  LogDebug("RecoMuon|CaloExtractorByAssociator")
314  << " >>> Veto HO tower: Calo deltaR= " << deltarHOcal;
315  LogDebug("RecoMuon|CaloExtractorByAssociator")
316  << " >>> Calo eta phi ethcal: " << calCPtr->eta() << " " << calCPtr->phi() << " " << ethcal;
317  vetoTowerHOCal = true;
318  }
319 
321  if (! (vetoTowerHOCal && vetoTowerHcal && vetoTowerEcal )){
322  for (unsigned int iH = 0; iH< mInfo.crossedTowerIds.size(); ++iH){
323  if (mInfo.crossedTowerIds[iH].rawId() == calCPtr->id().rawId()){
324  vetoTowerEcal = true;
325  vetoTowerHcal = true;
326  vetoTowerHOCal = true;
327  break;
328  }
329  }
330  }
331 
332  reco::isodeposit::Direction towerDir(calCPtr->eta(), calCPtr->phi());
334  if (doEcal){
335  if (vetoTowerEcal) depEcal.addCandEnergy(etecal);
336  else depEcal.addDeposit(towerDir, etecal);
337  }
338  if (doHcal){
339  if (vetoTowerHcal) depHcal.addCandEnergy(ethcal);
340  else depHcal.addDeposit(towerDir, ethcal);
341  }
342  if (doHOcal){
343  if (vetoTowerHOCal) depHOcal.addCandEnergy(ethocal);
344  else depHOcal.addDeposit(towerDir, ethocal);
345  }
346  }
347  }
348 
349  std::vector<IsoDeposit> resultDeps;
350  resultDeps.push_back(depEcal);
351  resultDeps.push_back(depHcal);
352  resultDeps.push_back(depHOcal);
353 
354  return resultDeps;
355 
356 }
357 
359  double noise = theNoiseTow_EB;
360  double eta = tower.eta();
361  if (fabs(eta)>1.479) noise = theNoiseTow_EE;
362  return noise;
363 }
364 
366  double noise = fabs(tower.eta())> 1.479 ? theNoise_HE : theNoise_HB;
367  return noise;
368 }
369 
371  double noise = theNoise_HO;
372  return noise;
373 }
374 
375 
376 double CaloExtractorByAssociator::noiseRecHit(const DetId& detId) const {
377  double noise = 100;
378  DetId::Detector det = detId.det();
379  if (det == DetId::Ecal){
380  EcalSubdetector subDet = (EcalSubdetector)(detId.subdetId());
381  if (subDet == EcalBarrel){
382  noise = theNoise_EB;
383  } else if (subDet == EcalEndcap){
384  noise = theNoise_EE;
385  }
386  } else if (det == DetId::Hcal){
387  HcalSubdetector subDet = (HcalSubdetector)(detId.subdetId());
388  if (subDet == HcalBarrel){
389  noise = theNoise_HB;
390  } else if (subDet == HcalEndcap){
391  noise = theNoise_HE;
392  } else if (subDet == HcalOuter){
393  noise = theNoise_HO;
394  }
395  }
396  return noise;
397 }
398 
#define LogDebug(id)
T getParameter(std::string const &) const
virtual void fillVetos(const edm::Event &ev, const edm::EventSetup &evSetup, const reco::TrackCollection &tracks)
allows to set extra vetoes (in addition to the muon) – no-op at this point
math::XYZPoint trkGlobPosAtHO
const Veto & veto() const
Get veto area.
Definition: IsoDeposit.h:73
std::vector< const CaloTower * > towers
virtual std::vector< reco::IsoDeposit > deposits(const edm::Event &ev, const edm::EventSetup &evSetup, const reco::Track &track) const
return deposits for 3 calorimeter subdetectors (ecal, hcal, ho) – in this order
void setVeto(const Veto &aVeto)
Set veto.
Definition: IsoDeposit.h:75
std::vector< DetId > crossedTowerIds
std::string thePropagatorName
propagator name to feed into the track associator
const DetId & detid() const
Definition: CaloRecHit.h:22
double hadEt() const
Definition: CaloTower.h:85
std::vector< const EcalRecHit * > ecalRecHits
hits in the cone
std::vector< DetId > crossedEcalIds
double theDR_Max
max cone size in which towers are considered
double outerEt() const
Definition: CaloTower.h:86
Geom::Phi< T > phi() const
Definition: PV3DBase.h:68
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:10
double noiseHOcal(const CaloTower &tower) const
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:139
std::vector< std::string > theDepositInstanceLabels
multiple deposits: labels – expect 3 labels beginning with &quot;e&quot;, &quot;h&quot;, &quot;ho&quot;
double noiseHcal(const CaloTower &tower) const
double theDR_Veto_E
cone sizes inside which the Et (towers) are not counted
TrackAssociatorParameters * theAssociatorParameters
associator, its&#39; parameters and the propagator
T eta() const
std::vector< DetId > crossedHcalIds
virtual double eta() const
momentum pseudorapidity
void addDeposit(double dr, double deposit)
Add deposit (ie. transverse energy or pT)
Definition: IsoDeposit.cc:23
math::XYZPoint trkGlobPosAtHcal
double deltaR(double eta1, double phi1, double eta2, double phi2)
Definition: deltaR.h:19
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
tuple TrackAssociatorParameters
Definition: example_cfg.py:32
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:141
double noiseEcal(const CaloTower &tower) const
Determine noise for HCAL and ECAL (take some defaults for the time being)
double emEnergy() const
Definition: CaloTower.h:79
void setPropagator(const Propagator *)
use a user configured propagator
float energy() const
Definition: CaloRecHit.h:19
std::vector< const HBHERecHit * > hcalRecHits
void addCandEnergy(double et)
Set energy or pT attached to cand trajectory.
Definition: IsoDeposit.h:134
std::vector< DetId > crossedHOIds
HcalSubdetector
Definition: HcalAssistant.h:32
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:39
double noiseRecHit(const DetId &detId) const
double hadEnergy() const
Definition: CaloTower.h:80
FreeTrajectoryState initialFreeState() const
Definition: DetId.h:20
CaloTowerDetId id() const
Definition: CaloTower.h:72
reco::IsoDeposit IsoDeposit
Definition: Isolation.h:7
static TimingReport * current()
Definition: TimingReport.cc:21
MuonServiceProxy * theService
the event setup proxy, it takes care the services update
Detector
Definition: DetId.h:26
const T & get() const
Definition: EventSetup.h:55
T eta() const
Definition: PV3DBase.h:75
tuple muons
Definition: patZpeak.py:38
TrackDetMatchInfo associate(const edm::Event &, const edm::EventSetup &, const FreeTrajectoryState &, const AssociatorParameters &)
math::XYZPoint trkGlobPosAtEcal
Track position at different parts of the calorimeter.
tuple cout
Definition: gather_cfg.py:121
virtual reco::IsoDeposit deposit(const edm::Event &ev, const edm::EventSetup &evSetup, const reco::Track &track) const
no-op: by design of this extractor the deposits are pulled out all at a time
bool thePrintTimeReport
flag to turn on/off printing of a time report
std::vector< const HORecHit * > hoRecHits
Detector det() const
get the detector field from this detid
Definition: DetId.h:37
EcalSubdetector
virtual double phi() const
momentum azimuthal angle
double emEt() const
Definition: CaloTower.h:84
void dump(std::ostream &co, bool active=false)
Definition: TimingReport.cc:50
double outerEnergy() const
Definition: CaloTower.h:81