CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HLTmmkFilter.cc
Go to the documentation of this file.
1 #include <algorithm>
2 #include <cmath>
3 #include <vector>
4 
8 
11 
14 
23 
25 
27 
29 
30 using namespace edm;
31 using namespace reco;
32 using namespace std;
33 using namespace trigger;
34 
35 
36 // ----------------------------------------------------------------------
38  muCandTag_ (iConfig.getParameter<edm::InputTag>("MuCand")),
39  muCandToken_(consumes<reco::RecoChargedCandidateCollection>(muCandTag_)),
40  trkCandTag_ (iConfig.getParameter<edm::InputTag>("TrackCand")),
41  trkCandToken_(consumes<reco::RecoChargedCandidateCollection>(trkCandTag_)),
42  thirdTrackMass_(iConfig.getParameter<double>("ThirdTrackMass")),
43  maxEta_(iConfig.getParameter<double>("MaxEta")),
44  minPt_(iConfig.getParameter<double>("MinPt")),
45  minInvMass_(iConfig.getParameter<double>("MinInvMass")),
46  maxInvMass_(iConfig.getParameter<double>("MaxInvMass")),
47  maxNormalisedChi2_(iConfig.getParameter<double>("MaxNormalisedChi2")),
48  minLxySignificance_(iConfig.getParameter<double>("MinLxySignificance")),
49  minCosinePointingAngle_(iConfig.getParameter<double>("MinCosinePointingAngle")),
50  minD0Significance_(iConfig.getParameter<double>("MinD0Significance")),
51  fastAccept_(iConfig.getParameter<bool>("FastAccept")),
52  beamSpotTag_ (iConfig.getParameter<edm::InputTag> ("BeamSpotTag")),
53  beamSpotToken_(consumes<reco::BeamSpot>(beamSpotTag_))
54 {
55  produces<VertexCollection>();
56  produces<CandidateCollection>();
57 }
58 
59 
60 // ----------------------------------------------------------------------
62 
63 }
64 
65 void
69  desc.add<edm::InputTag>("MuCand",edm::InputTag("hltMuTracks"));
70  desc.add<edm::InputTag>("TrackCand",edm::InputTag("hltMumukAllConeTracks"));
71  desc.add<double>("ThirdTrackMass",0.106);
72  desc.add<double>("MaxEta",2.5);
73  desc.add<double>("MinPt",3.0);
74  desc.add<double>("MinInvMass",1.2);
75  desc.add<double>("MaxInvMass",2.2);
76  desc.add<double>("MaxNormalisedChi2",10.0);
77  desc.add<double>("MinLxySignificance",3.0);
78  desc.add<double>("MinCosinePointingAngle",0.9);
79  desc.add<double>("MinD0Significance",0.0);
80  desc.add<bool>("FastAccept",false);
81  desc.add<edm::InputTag>("BeamSpotTag",edm::InputTag("hltOfflineBeamSpot"));
82  descriptions.add("hltmmkFilter",desc);
83 }
84 
85 // ----------------------------------------------------------------------
87 
88 }
89 
90 
91 // ----------------------------------------------------------------------
93 
94 }
95 
96 
97 // ----------------------------------------------------------------------
99 
100  const double MuMass(0.106);
101  const double MuMass2(MuMass*MuMass);
102 
103  const double thirdTrackMass2(thirdTrackMass_*thirdTrackMass_);
104 
105  auto_ptr<CandidateCollection> output(new CandidateCollection());
106  auto_ptr<VertexCollection> vertexCollection(new VertexCollection());
107 
108  //get the transient track builder:
110  iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theB);
111 
112  //get the beamspot position
113  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
114  iEvent.getByToken(beamSpotToken_,recoBeamSpotHandle);
115  const reco::BeamSpot& vertexBeamSpot = *recoBeamSpotHandle;
116 
117  ESHandle<MagneticField> bFieldHandle;
118  iSetup.get<IdealMagneticFieldRecord>().get(bFieldHandle);
119 
120  const MagneticField* magField = bFieldHandle.product();
121 
122  TSCBLBuilderNoMaterial blsBuilder;
123 
124  // Ref to Candidate object to be recorded in filter object
127  RecoChargedCandidateRef refTrk;
128 
129  // get hold of muon trks
131  iEvent.getByToken(muCandToken_,mucands);
132 
133  // get track candidates around displaced muons
135  iEvent.getByToken(trkCandToken_,trkcands);
136 
137  if (saveTags()) {
138  filterproduct.addCollectionTag(muCandTag_);
139  filterproduct.addCollectionTag(trkCandTag_);
140  }
141 
142  double e1,e2,e3;
144 
145  //TrackRefs to mu cands in trkcand
146  vector<TrackRef> trkMuCands;
147 
148  //Already used mu tracks to avoid double counting candidates
149  vector<bool> isUsedCand(trkcands->size(),false);
150 
151  int counter = 0;
152 
153  //run algorithm
154  for (RecoChargedCandidateCollection::const_iterator mucand1=mucands->begin(), endCand1=mucands->end(); mucand1!=endCand1; ++mucand1) {
155 
156  if ( mucands->size()<2) break;
157  if ( trkcands->size()<1) break;
158 
159  TrackRef trk1 = mucand1->get<TrackRef>();
160  LogDebug("HLTDisplacedMumukFilter") << " 1st muon: q*pt= " << trk1->charge()*trk1->pt() << ", eta= " << trk1->eta() << ", hits= " << trk1->numberOfValidHits();
161 
162  // eta cut
163  if (fabs(trk1->eta()) > maxEta_) continue;
164 
165  // Pt threshold cut
166  if (trk1->pt() < minPt_) continue;
167 
168  RecoChargedCandidateCollection::const_iterator mucand2 = mucand1; ++mucand2;
169 
170  for (RecoChargedCandidateCollection::const_iterator endCand2=mucands->end(); mucand2!=endCand2; ++mucand2) {
171 
172  TrackRef trk2 = mucand2->get<TrackRef>();
173 
174  LogDebug("HLTDisplacedMumukFilter") << " 2nd muon: q*pt= " << trk2->charge()*trk2->pt() << ", eta= " << trk2->eta() << ", hits= " << trk2->numberOfValidHits();
175 
176  // eta cut
177  if (fabs(trk2->eta()) > maxEta_) continue;
178 
179  // Pt threshold cut
180  if (trk2->pt() < minPt_) continue;
181 
182  RecoChargedCandidateCollection::const_iterator trkcand, endCandTrk;
183 
184  std::vector<bool>::iterator isUsedIter, endIsUsedCand;
185 
186  //get overlapping muon candidates
187  for ( trkcand = trkcands->begin(), endCandTrk=trkcands->end(), isUsedIter = isUsedCand.begin(), endIsUsedCand = isUsedCand.end(); trkcand != endCandTrk && isUsedIter != endIsUsedCand; ++trkcand, ++isUsedIter) {
188  TrackRef trk3 = trkcand->get<TrackRef>();
189 
190  //check for overlapping muon tracks and save TrackRefs
191  if (overlap(*mucand1,*trkcand)) {
192  trkMuCands.push_back(trk3);
193  *isUsedIter = true;
194  continue;
195  }
196  else if (overlap(*mucand2,*trkcand)){
197  trkMuCands.push_back(trk3);
198  *isUsedIter = true;
199  continue;
200  }
201 
202  if(trkMuCands.size()==2) break;
203  }
204 
205  //Not all overlapping candidates found, skip event
206  //if (trkMuCands.size()!=2) continue;
207 
208  //combine muons with all tracks
209  for ( trkcand = trkcands->begin(), endCandTrk=trkcands->end(), isUsedIter = isUsedCand.begin(), endIsUsedCand = isUsedCand.end(); trkcand != endCandTrk && isUsedIter != endIsUsedCand; ++trkcand, ++isUsedIter) {
210 
211  TrackRef trk3 = trkcand->get<TrackRef>();
212 
213  LogDebug("HLTDisplacedMumukFilter") << " 3rd track: q*pt= " << trk3->charge()*trk3->pt() << ", eta= " << trk3->eta() << ", hits= " << trk3->numberOfValidHits();
214 
215  //skip overlapping muon candidates
216  bool skip=false;
217  for (unsigned int itmc=0;itmc<trkMuCands.size();itmc++) if(trk3==trkMuCands.at(itmc)) skip=true;
218  if(skip) continue;
219 
220  //skip already used tracks
221  if(*isUsedIter) continue;
222 
223  // eta cut
224  if (fabs(trk3->eta()) > maxEta_) continue;
225 
226  // Pt threshold cut
227  if (trk3->pt() < minPt_) continue;
228 
229  // Combined system
230  e1 = sqrt(trk1->momentum().Mag2()+MuMass2);
231  e2 = sqrt(trk2->momentum().Mag2()+MuMass2);
232  e3 = sqrt(trk3->momentum().Mag2()+thirdTrackMass2);
233 
234  p1 = Particle::LorentzVector(trk1->px(),trk1->py(),trk1->pz(),e1);
235  p2 = Particle::LorentzVector(trk2->px(),trk2->py(),trk2->pz(),e2);
236  p3 = Particle::LorentzVector(trk3->px(),trk3->py(),trk3->pz(),e3);
237 
238  p = p1+p2+p3;
239 
240  //invariant mass cut
241  double invmass = abs(p.mass());
242 
243  LogDebug("HLTDisplacedMumukFilter") << " Invmass= " << invmass;
244 
245  if (invmass>maxInvMass_ || invmass<minInvMass_) continue;
246 
247  // do the vertex fit
248  vector<TransientTrack> t_tks;
249  t_tks.push_back((*theB).build(&trk1));
250  t_tks.push_back((*theB).build(&trk2));
251  t_tks.push_back((*theB).build(&trk3));
252 
253  if (t_tks.size()!=3) continue;
254 
255  FreeTrajectoryState InitialFTS = initialFreeState(*trk3, magField);
256  TrajectoryStateClosestToBeamLine tscb( blsBuilder(InitialFTS, *recoBeamSpotHandle) );
257  double d0sig = tscb.transverseImpactParameter().significance();
258 
259  if (d0sig < minD0Significance_) continue;
260 
261  KalmanVertexFitter kvf;
262  TransientVertex tv = kvf.vertex(t_tks);
263 
264  if (!tv.isValid()) continue;
265 
266  Vertex vertex = tv;
267 
268  // get vertex position and error to calculate the decay length significance
269  GlobalPoint secondaryVertex = tv.position();
270  GlobalError err = tv.positionError();
271 
272  //calculate decay length significance w.r.t. the beamspot
273  GlobalPoint displacementFromBeamspot( -1*((vertexBeamSpot.x0() -secondaryVertex.x()) + (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dxdz()), -1*((vertexBeamSpot.y0() - secondaryVertex.y())+ (secondaryVertex.z() -vertexBeamSpot.z0()) * vertexBeamSpot.dydz()), 0);
274 
275  float lxy = displacementFromBeamspot.perp();
276  float lxyerr = sqrt(err.rerr(displacementFromBeamspot));
277 
278  // get normalizes chi2
279  float normChi2 = tv.normalisedChiSquared();
280 
281  //calculate the angle between the decay length and the mumu momentum
282  Vertex::Point vperp(displacementFromBeamspot.x(),displacementFromBeamspot.y(),0.);
283  math::XYZVector pperp(p.x(),p.y(),0.);
284 
285  float cosAlpha = vperp.Dot(pperp)/(vperp.R()*pperp.R());
286 
287  LogDebug("HLTDisplacedMumukFilter") << " vertex fit normalised chi2: " << normChi2 << ", Lxy significance: " << lxy/lxyerr << ", cosine pointing angle: " << cosAlpha;
288 
289  // put vertex in the event
290  vertexCollection->push_back(vertex);
291 
292  if (normChi2 > maxNormalisedChi2_) continue;
293  if (lxy/lxyerr < minLxySignificance_) continue;
294  if(cosAlpha < minCosinePointingAngle_) continue;
295 
296  LogDebug("HLTDisplacedMumukFilter") << " Event passed!";
297 
298  //Add event
299  ++counter;
300 
301  //Get refs
302  bool i1done = false;
303  bool i2done = false;
304  bool i3done = false;
305  vector<RecoChargedCandidateRef> vref;
306  filterproduct.getObjects(TriggerMuon,vref);
307  for (unsigned int i=0; i<vref.size(); i++) {
309  TrackRef trktmp = candref->get<TrackRef>();
310  if (trktmp==trk1) {
311  i1done = true;
312  } else if (trktmp==trk2) {
313  i2done = true;
314  } else if (trktmp==trk3) {
315  i3done = true;
316  }
317  if (i1done && i2done && i3done) break;
318  }
319 
320  if (!i1done) {
321  refMu1=RecoChargedCandidateRef( Ref<RecoChargedCandidateCollection> (mucands,distance(mucands->begin(), mucand1)));
322  filterproduct.addObject(TriggerMuon,refMu1);
323  }
324  if (!i2done) {
325  refMu2=RecoChargedCandidateRef( Ref<RecoChargedCandidateCollection> (mucands,distance(mucands->begin(),mucand2)));
326  filterproduct.addObject(TriggerMuon,refMu2);
327  }
328  if (!i3done) {
329  refTrk=RecoChargedCandidateRef( Ref<RecoChargedCandidateCollection> (trkcands,distance(trkcands->begin(),trkcand)));
330  filterproduct.addObject(TriggerTrack,refTrk);
331  }
332 
333  if (fastAccept_) break;
334  }
335 
336  trkMuCands.clear();
337  }
338  }
339 
340  // filter decision
341  const bool accept (counter >= 1);
342 
343  LogDebug("HLTDisplacedMumukFilter") << " >>>>> Result of HLTDisplacedMumukFilter is "<< accept << ", number of muon pairs passing thresholds= " << counter;
344 
345  iEvent.put(vertexCollection);
346 
347  return accept;
348 }
349 
350 // ----------------------------------------------------------------------
352  const MagneticField* field)
353 {
354  Basic3DVector<float> pos( tk.vertex());
355  GlobalPoint gpos( pos);
356  Basic3DVector<float> mom( tk.momentum());
357  GlobalVector gmom( mom);
358  GlobalTrajectoryParameters par( gpos, gmom, tk.charge(), field);
360  return FreeTrajectoryState( par, err);
361 }
362 
364 
365  double eps(1.44e-4);
366 
367  double dpt = a.pt() - b.pt();
368  dpt *= dpt;
369 
370  double dphi = deltaPhi(a.phi(), b.phi());
371  dphi *= dphi;
372 
373  double deta = a.eta() - b.eta();
374  deta *= deta;
375 
376  if ((dphi + deta) < eps) {
377  return 1;
378  }
379 
380  return 0;
381 
382 }
383 
#define LogDebug(id)
GlobalError positionError() const
static int overlap(const reco::Candidate &, const reco::Candidate &)
double z0() const
z coordinate
Definition: BeamSpot.h:68
int i
Definition: DBlmapReader.cc:9
edm::EDGetTokenT< reco::RecoChargedCandidateCollection > trkCandToken_
Definition: HLTmmkFilter.h:59
HLTmmkFilter(const edm::ParameterSet &)
Definition: HLTmmkFilter.cc:37
void getObjects(Vids &ids, VRphoton &refs) const
various physics-level getters:
T perp() const
Definition: PV3DBase.h:72
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: HLTmmkFilter.cc:66
const double maxInvMass_
Definition: HLTmmkFilter.h:65
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:464
virtual double pt() const =0
transverse momentum
virtual CachingVertex< 5 > vertex(const std::vector< reco::TransientTrack > &tracks) const
T y() const
Definition: PV3DBase.h:63
edm::Ref< RecoChargedCandidateCollection > RecoChargedCandidateRef
reference to an object in a collection of RecoChargedCandidate objects
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
const Vector & momentum() const
track momentum vector
Definition: TrackBase.h:662
bool accept(const edm::Event &event, const edm::TriggerResults &triggerTable, const std::string &triggerPath)
Definition: TopDQMHelpers.h:24
tuple vertexCollection
virtual bool hltFilter(edm::Event &, const edm::EventSetup &, trigger::TriggerFilterObjectWithRefs &filterproduct) const override
Definition: HLTmmkFilter.cc:98
void addObject(int id, const reco::RecoEcalCandidateRef &ref)
setters for L3 collections: (id=physics type, and Ref&lt;C&gt;)
double dydz() const
dydz slope
Definition: BeamSpot.h:84
float normalisedChiSquared() const
const Point & vertex() const
reference point on the track. This method is DEPRECATED, please use referencePoint() instead ...
Definition: TrackBase.h:674
int iEvent
Definition: GenABIO.cc:230
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:120
CovarianceMatrix covariance() const
return track covariance matrix
Definition: TrackBase.h:718
T sqrt(T t)
Definition: SSEVec.h:48
GlobalPoint position() const
T z() const
Definition: PV3DBase.h:64
virtual void beginJob()
Definition: HLTmmkFilter.cc:86
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
math::XYZPoint Point
point in the space
Definition: Vertex.h:39
edm::EDGetTokenT< reco::BeamSpot > beamSpotToken_
Definition: HLTmmkFilter.h:72
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:244
edm::InputTag trkCandTag_
Definition: HLTmmkFilter.h:58
ParameterDescriptionBase * add(U const &iLabel, T const &value)
double p2[4]
Definition: TauolaWrapper.h:90
double dxdz() const
dxdz slope
Definition: BeamSpot.h:82
edm::EDGetTokenT< reco::RecoChargedCandidateCollection > muCandToken_
Definition: HLTmmkFilter.h:57
std::vector< RecoChargedCandidate > RecoChargedCandidateCollection
collectin of RecoChargedCandidate objects
double significance() const
Definition: Measurement1D.h:32
static void makeHLTFilterDescription(edm::ParameterSetDescription &desc)
Definition: HLTFilter.cc:29
edm::InputTag muCandTag_
Definition: HLTmmkFilter.h:56
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
const double maxEta_
Definition: HLTmmkFilter.h:62
T rerr(const GlobalPoint &aPoint) const
const T & get() const
Definition: EventSetup.h:56
T const * product() const
Definition: ESHandle.h:86
double b
Definition: hdecay.h:120
void addCollectionTag(const edm::InputTag &collectionTag)
collectionTags
void add(std::string const &label, ParameterSetDescription const &psetDescription)
virtual void endJob()
Definition: HLTmmkFilter.cc:92
const double minCosinePointingAngle_
Definition: HLTmmkFilter.h:68
const double minInvMass_
Definition: HLTmmkFilter.h:64
const double thirdTrackMass_
Definition: HLTmmkFilter.h:61
bool saveTags() const
Definition: HLTFilter.h:45
double p1[4]
Definition: TauolaWrapper.h:89
double a
Definition: hdecay.h:121
static std::atomic< unsigned int > counter
double y0() const
y coordinate
Definition: BeamSpot.h:66
int charge() const
track electric charge
Definition: TrackBase.h:554
const double minPt_
Definition: HLTmmkFilter.h:63
const double maxNormalisedChi2_
Definition: HLTmmkFilter.h:66
static FreeTrajectoryState initialFreeState(const reco::Track &, const MagneticField *)
edm::OwnVector< Candidate > CandidateCollection
collection of Candidate objects
Definition: CandidateFwd.h:21
T x() const
Definition: PV3DBase.h:62
const double minD0Significance_
Definition: HLTmmkFilter.h:69
bool isValid() const
math::PtEtaPhiELorentzVectorF LorentzVector
const bool fastAccept_
Definition: HLTmmkFilter.h:70
virtual double phi() const =0
momentum azimuthal angle
virtual double eta() const =0
momentum pseudorapidity
const double minLxySignificance_
Definition: HLTmmkFilter.h:67
double p3[4]
Definition: TauolaWrapper.h:91
double x0() const
x coordinate
Definition: BeamSpot.h:64