test
CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TauDiscriminationAgainstCaloMuon.cc
Go to the documentation of this file.
1 
2 /*
3  * class TauDiscriminationAgainstCaloMuon
4  * created : Nov 20 2010
5  * revised :
6  * Authors : Christian Veelken (UC Davis)
7  */
8 
10 
13 
16 
24 
29 
34 
37 
39 
40 #include <TVector3.h>
41 #include <TMath.h>
42 
43 #include <string>
44 
45 namespace {
46 
47 using namespace reco;
48 
49 // define acess to lead. track for CaloTaus
50 template <typename T>
51 class TauLeadTrackExtractor
52 {
53  public:
54  reco::TrackRef getLeadTrack(const T& tau) const
55  {
56  return tau.leadTrack();
57  }
58  double getTrackPtSum(const T& tau) const
59  {
60  double trackPtSum = 0.;
61  for ( TrackRefVector::const_iterator signalTrack = tau.signalTracks().begin();
62  signalTrack != tau.signalTracks().end(); ++signalTrack ) {
63  trackPtSum += (*signalTrack)->pt();
64  }
65  return trackPtSum;
66  }
67 };
68 
69 // define acess to lead. track for PFTaus
70 template <>
71 class TauLeadTrackExtractor<reco::PFTau>
72 {
73  public:
74  reco::TrackRef getLeadTrack(const reco::PFTau& tau) const
75  {
76  return tau.leadPFChargedHadrCand()->trackRef();
77  }
78  double getTrackPtSum(const reco::PFTau& tau) const
79  {
80  double trackPtSum = 0.;
81  for ( std::vector<PFCandidatePtr>::const_iterator signalTrack = tau.signalPFChargedHadrCands().begin();
82  signalTrack != tau.signalPFChargedHadrCands().end(); ++signalTrack ) {
83  trackPtSum += (*signalTrack)->pt();
84  }
85  return trackPtSum;
86  }
87 };
88 
89 template<class TauType, class TauDiscriminator>
90 class TauDiscriminationAgainstCaloMuon final : public TauDiscriminationProducerBase<TauType, TauDiscriminator>
91 {
92  public:
93  // setup framework types for this tautype
94  typedef std::vector<TauType> TauCollection;
96 
97  explicit TauDiscriminationAgainstCaloMuon(const edm::ParameterSet&);
98  ~TauDiscriminationAgainstCaloMuon() {}
99 
100  // called at the beginning of every event
101  void beginEvent(const edm::Event&, const edm::EventSetup&);
102 
103  double discriminate(const TauRef&) const override;
104 
105  private:
106  edm::InputTag srcEcalRecHitsBarrel_;
108  edm::InputTag srcEcalRecHitsEndcap_;
110  edm::InputTag srcHcalRecHits_;
112 
113  edm::InputTag srcVertex_;
114  GlobalPoint eventVertexPosition_;
115 
116  const TransientTrackBuilder* trackBuilder_;
117  const CaloGeometry* caloGeometry_;
118 
119  TauLeadTrackExtractor<TauType> leadTrackExtractor_;
120 
121  double minLeadTrackPt_;
122  double minLeadTrackPtFraction_;
123 
124  double drEcal_;
125  double drHcal_;
126 
127  double maxEnEcal_;
128  double maxEnHcal_;
129 
130  double maxEnToTrackRatio_;
131 };
132 
133 template<class TauType, class TauDiscriminator>
134 TauDiscriminationAgainstCaloMuon<TauType, TauDiscriminator>::TauDiscriminationAgainstCaloMuon(const edm::ParameterSet& cfg)
135  : TauDiscriminationProducerBase<TauType, TauDiscriminator>(cfg)
136 {
137  srcEcalRecHitsBarrel_ = cfg.getParameter<edm::InputTag>("srcEcalRecHitsBarrel");
138  srcEcalRecHitsEndcap_ = cfg.getParameter<edm::InputTag>("srcEcalRecHitsEndcap");
139  srcHcalRecHits_ = cfg.getParameter<edm::InputTag>("srcHcalRecHits");
140 
141  srcVertex_ = cfg.getParameter<edm::InputTag>("srcVertex");
142 
143  minLeadTrackPt_ = cfg.getParameter<double>("minLeadTrackPt");
144  minLeadTrackPtFraction_ = cfg.getParameter<double>("minLeadTrackPtFraction");
145 
146  drEcal_ = cfg.getParameter<double>("dRecal");
147  drHcal_ = cfg.getParameter<double>("dRhcal");
148 
149  maxEnEcal_ = cfg.getParameter<double>("maxEnEcal");
150  maxEnHcal_ = cfg.getParameter<double>("maxEnHcal");
151 
152  maxEnToTrackRatio_ = cfg.getParameter<double>("maxEnToTrackRatio");
153 }
154 
155 template<class TauType, class TauDiscriminator>
156 void TauDiscriminationAgainstCaloMuon<TauType, TauDiscriminator>::beginEvent(const edm::Event& evt, const edm::EventSetup& evtSetup)
157 {
158  evt.getByLabel(srcEcalRecHitsBarrel_, ebRecHits_);
159  evt.getByLabel(srcEcalRecHitsEndcap_, eeRecHits_);
160  evt.getByLabel(srcHcalRecHits_, hbheRecHits_);
161 
162  edm::ESHandle<TransientTrackBuilder> trackBuilderHandle;
163  evtSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", trackBuilderHandle);
164  trackBuilder_ = trackBuilderHandle.product();
165  if ( !trackBuilder_ ) {
166  edm::LogError ("TauDiscriminationAgainstCaloMuon::discriminate")
167  << " Failed to access TransientTrackBuilder !!";
168  }
169 
170  edm::ESHandle<CaloGeometry> caloGeometryHandle;
171  evtSetup.get<CaloGeometryRecord>().get(caloGeometryHandle);
172  caloGeometry_ = caloGeometryHandle.product();
173  if ( !caloGeometry_ ) {
174  edm::LogError ("TauDiscriminationAgainstCaloMuon::discriminate")
175  << " Failed to access CaloGeometry !!";
176  }
177 
179  evt.getByLabel(srcVertex_, vertices);
180  eventVertexPosition_ = GlobalPoint(0., 0., 0.);
181  if ( vertices->size() >= 1 ) {
182  const reco::Vertex& thePrimaryEventVertex = (*vertices->begin());
183  eventVertexPosition_ = GlobalPoint(thePrimaryEventVertex.x(), thePrimaryEventVertex.y(), thePrimaryEventVertex.z());
184  }
185 }
186 
187 double compEcalEnergySum(const EcalRecHitCollection& ecalRecHits,
188  const CaloSubdetectorGeometry* detGeometry,
189  const reco::TransientTrack& transientTrack, double dR,
190  const GlobalPoint& eventVertexPosition)
191 {
192  double ecalEnergySum = 0.;
194  ecalRecHit != ecalRecHits.end(); ++ecalRecHit ) {
195  const CaloCellGeometry* cellGeometry = detGeometry->getGeometry(ecalRecHit->detid());
196 
197  if ( !cellGeometry ) {
198  edm::LogError ("compEcalEnergySum")
199  << " Failed to access ECAL geometry for detId = " << ecalRecHit->detid().rawId()
200  << " --> skipping !!";
201  continue;
202  }
203 
204  const GlobalPoint& cellPosition = cellGeometry->getPosition();
205 
206 //--- CV: speed up computation by requiring eta-phi distance
207 // between cell position and track direction to be dR < 0.5
208  Vector3DBase<float, GlobalTag> cellPositionRelVertex = (cellPosition) - eventVertexPosition;
209  if ( deltaR(cellPositionRelVertex.eta(), cellPositionRelVertex.phi(),
210  transientTrack.track().eta(), transientTrack.track().phi()) > 0.5 ) continue;
211 
212  TrajectoryStateClosestToPoint dcaPosition = transientTrack.trajectoryStateClosestToPoint(cellPosition);
213 
214  Vector3DBase<float, GlobalTag> d = (cellPosition - dcaPosition.position());
215 
216  TVector3 d3(d.x(), d.y(), d.z());
217  TVector3 dir(transientTrack.track().px(), transientTrack.track().py(), transientTrack.track().pz());
218 
219  double dPerp = d3.Cross(dir.Unit()).Mag();
220  double dParl = TVector3(cellPosition.x(), cellPosition.y(), cellPosition.z()).Dot(dir.Unit());
221 
222  if ( dPerp < dR && dParl > 100. ) {
223  ecalEnergySum += ecalRecHit->energy();
224  }
225  }
226 
227  return ecalEnergySum;
228 }
229 
230 double compHcalEnergySum(const HBHERecHitCollection& hcalRecHits,
231  const CaloSubdetectorGeometry* hbGeometry, const CaloSubdetectorGeometry* heGeometry,
232  const reco::TransientTrack& transientTrack, double dR,
233  const GlobalPoint& eventVertexPosition)
234 {
235  double hcalEnergySum = 0.;
236  for ( HBHERecHitCollection::const_iterator hcalRecHit = hcalRecHits.begin();
237  hcalRecHit != hcalRecHits.end(); ++hcalRecHit ) {
238  const CaloCellGeometry* hbCellGeometry = hbGeometry->getGeometry(hcalRecHit->detid());
239  const CaloCellGeometry* heCellGeometry = heGeometry->getGeometry(hcalRecHit->detid());
240 
241  const GlobalPoint* cellPosition = 0;
242  if ( hbCellGeometry ) cellPosition = &(hbCellGeometry->getPosition());
243  if ( heCellGeometry ) cellPosition = &(heCellGeometry->getPosition());
244 
245  if ( !cellPosition ) {
246  edm::LogError ("compHcalEnergySum")
247  << " Failed to access HCAL geometry for detId = " << hcalRecHit->detid().rawId()
248  << " --> skipping !!";
249  continue;
250  }
251 
252 //--- CV: speed up computation by requiring eta-phi distance
253 // between cell position and track direction to be dR < 0.5
254  Vector3DBase<float, GlobalTag> cellPositionRelVertex = (*cellPosition) - eventVertexPosition;
255  if ( deltaR(cellPositionRelVertex.eta(), cellPositionRelVertex.phi(),
256  transientTrack.track().eta(), transientTrack.track().phi()) > 0.5 ) continue;
257 
258  TrajectoryStateClosestToPoint dcaPosition = transientTrack.trajectoryStateClosestToPoint(*cellPosition);
259 
260  Vector3DBase<float, GlobalTag> d = ((*cellPosition) - dcaPosition.position());
261 
262  TVector3 d3(d.x(), d.y(), d.z());
263  TVector3 dir(transientTrack.track().px(), transientTrack.track().py(), transientTrack.track().pz());
264 
265  double dPerp = d3.Cross(dir.Unit()).Mag();
266  double dParl = TVector3(cellPosition->x(), cellPosition->y(), cellPosition->z()).Dot(dir.Unit());
267 
268  if ( dPerp < dR && dParl > 100. ) {
269  hcalEnergySum += hcalRecHit->energy();
270  }
271  }
272 
273  return hcalEnergySum;
274 }
275 
276 template<class TauType, class TauDiscriminator>
277 double TauDiscriminationAgainstCaloMuon<TauType, TauDiscriminator>::discriminate(const TauRef& tau) const
278 {
279  if ( !(trackBuilder_ && caloGeometry_) ) return 0.;
280 
281  const CaloSubdetectorGeometry* ebGeometry = caloGeometry_->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
282  const CaloSubdetectorGeometry* eeGeometry = caloGeometry_->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
283  const CaloSubdetectorGeometry* hbGeometry = caloGeometry_->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
284  const CaloSubdetectorGeometry* heGeometry = caloGeometry_->getSubdetectorGeometry(DetId::Hcal, HcalEndcap);
285 
286  TrackRef leadTrackRef = leadTrackExtractor_.getLeadTrack(*tau);
287 
288  if ( (leadTrackRef.isAvailable() || leadTrackRef.isTransient()) && leadTrackRef.isNonnull() ) {
289  double leadTrackPt = leadTrackRef->pt();
290  double trackPtSum = leadTrackExtractor_.getTrackPtSum(*tau);
291 
292  double leadTrackPtFraction = ( trackPtSum > 0. ) ? (leadTrackPt/trackPtSum) : -1.;
293 
294  if ( leadTrackPt > minLeadTrackPt_ && leadTrackPtFraction > minLeadTrackPtFraction_ ) {
295  reco::TransientTrack transientTrack = trackBuilder_->build(leadTrackRef);
296 
297  double ebEnergySum = compEcalEnergySum(*ebRecHits_, ebGeometry, transientTrack, drEcal_, eventVertexPosition_);
298  double eeEnergySum = compEcalEnergySum(*eeRecHits_, eeGeometry, transientTrack, drEcal_, eventVertexPosition_);
299  double ecalEnergySum = ebEnergySum + eeEnergySum;
300 
301  double hbheEnergySum = compHcalEnergySum(*hbheRecHits_, hbGeometry, heGeometry, transientTrack, drHcal_, eventVertexPosition_);
302 
303  double caloEnergySum = ecalEnergySum + hbheEnergySum;
304 
305  if ( ecalEnergySum < maxEnEcal_ &&
306  hbheEnergySum < maxEnHcal_ &&
307  caloEnergySum < (maxEnToTrackRatio_*leadTrackPt) ) return 0.;
308  }
309  }
310 
311  return 1.;
312 }
313 
314 }
315 
317 
318 typedef TauDiscriminationAgainstCaloMuon<PFTau, PFTauDiscriminator> PFRecoTauDiscriminationAgainstCaloMuon;
319 typedef TauDiscriminationAgainstCaloMuon<CaloTau, CaloTauDiscriminator> CaloRecoTauDiscriminationAgainstCaloMuon;
320 
321 DEFINE_FWK_MODULE(PFRecoTauDiscriminationAgainstCaloMuon);
322 DEFINE_FWK_MODULE(CaloRecoTauDiscriminationAgainstCaloMuon);
bool isAvailable() const
Definition: Ref.h:576
T getParameter(std::string const &) const
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:252
tuple cfg
Definition: looper.py:293
edm::Ref< TauBxCollection > TauRef
Definition: Tau.h:12
const PFCandidatePtr & leadPFChargedHadrCand() const
Definition: PFTau.cc:67
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
double y() const
y coordinate
Definition: Vertex.h:103
bool isTransient() const
Checks if this ref is transient (i.e. not persistable).
Definition: Ref.h:277
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
std::vector< EcalRecHit >::const_iterator const_iterator
T y() const
Definition: PV3DBase.h:63
std::vector< Tau > TauCollection
Definition: Tau.h:38
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:640
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:622
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
tuple d
Definition: ztail.py:151
TauDiscriminationAgainstCaloMuon< PFTau, PFTauDiscriminator > PFRecoTauDiscriminationAgainstCaloMuon
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:646
T z() const
Definition: PV3DBase.h:64
double z() const
y coordinate
Definition: Vertex.h:105
auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:28
TauDiscriminationAgainstCaloMuon< CaloTau, CaloTauDiscriminator > CaloRecoTauDiscriminationAgainstCaloMuon
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:418
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:634
const_iterator end() const
double x() const
x coordinate
Definition: Vertex.h:101
const Track & track() const
TrajectoryStateClosestToPoint trajectoryStateClosestToPoint(const GlobalPoint &point) const
const T & get() const
Definition: EventSetup.h:56
T const * product() const
Definition: ESHandle.h:86
T eta() const
Definition: PV3DBase.h:76
dbl *** dir
Definition: mlp_gen.cc:35
long double T
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
T x() const
Definition: PV3DBase.h:62
const std::vector< reco::PFCandidatePtr > & signalPFChargedHadrCands() const
Charged hadrons in signal region.
Definition: PFTau.cc:80
const_iterator begin() const
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:628