CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CosmicRegionalSeedGenerator.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <memory>
3 #include <string>
4 
9 
13 
16 
20 
23 
26 
29 
31 
32 #include <TMath.h>
33 
34 using namespace std;
35 using namespace trigger;
36 using namespace reco;
37 using namespace edm;
38 
40  conf_(conf)
41 {
42  edm::LogInfo ("CosmicRegionalSeedGenerator") << "Begin Run:: Constructing CosmicRegionalSeedGenerator";
43 
45  ptMin_ = regionPSet.getParameter<double>("ptMin");
46  rVertex_ = regionPSet.getParameter<double>("rVertex");
47  zVertex_ = regionPSet.getParameter<double>("zVertex");
48  deltaEta_ = regionPSet.getParameter<double>("deltaEtaRegion");
49  deltaPhi_ = regionPSet.getParameter<double>("deltaPhiRegion");
50 
51  edm::ParameterSet toolsPSet = conf_.getParameter<edm::ParameterSet>("ToolsPSet");
52  thePropagatorName_ = toolsPSet.getParameter<std::string>("thePropagatorName");
53  regionBase_ = toolsPSet.getParameter<std::string>("regionBase");
54 
55  edm::ParameterSet collectionsPSet = conf_.getParameter<edm::ParameterSet>("CollectionsPSet");
56  recoMuonsCollection_ = collectionsPSet.getParameter<edm::InputTag>("recoMuonsCollection");
57  recoTrackMuonsCollection_ = collectionsPSet.getParameter<edm::InputTag>("recoTrackMuonsCollection");
58  recoL2MuonsCollection_ = collectionsPSet.getParameter<edm::InputTag>("recoL2MuonsCollection");
59 
60  edm::ParameterSet regionInJetsCheckPSet = conf_.getParameter<edm::ParameterSet>("RegionInJetsCheckPSet");
61  doJetsExclusionCheck_ = regionInJetsCheckPSet.getParameter<bool>("doJetsExclusionCheck");
62  deltaRExclusionSize_ = regionInJetsCheckPSet.getParameter<double>("deltaRExclusionSize");
63  jetsPtMin_ = regionInJetsCheckPSet.getParameter<double>("jetsPtMin");
64  recoCaloJetsCollection_ = regionInJetsCheckPSet.getParameter<edm::InputTag>("recoCaloJetsCollection");
65 
66 
67  edm::LogInfo ("CosmicRegionalSeedGenerator") << "Reco muons collection: " << recoMuonsCollection_ << "\n"
68  << "Reco tracks muons collection: " << recoTrackMuonsCollection_<< "\n"
69  << "Reco L2 muons collection: " << recoL2MuonsCollection_;
70 }
71 
72 std::vector<TrackingRegion*, std::allocator<TrackingRegion*> > CosmicRegionalSeedGenerator::regions(const edm::Event& event, const edm::EventSetup& es) const
73 {
74 
75  std::vector<TrackingRegion* > result;
76 
77 
78  //________________________________________
79  //
80  //Seeding on Sta muon (MC && Datas)
81  //________________________________________
82 
83 
84  if(regionBase_=="seedOnStaMuon"||regionBase_=="") {
85 
86  LogDebug("CosmicRegionalSeedGenerator") << "Seeding on stand alone muons ";
87 
88  //get collections
89  //+++++++++++++++
90 
91  //get the muon collection
93  event.getByLabel(recoMuonsCollection_,muonsHandle);
94  if (!muonsHandle.isValid())
95  {
96  edm::LogError("CollectionNotFound") << "Error::No reco muons collection (" << recoMuonsCollection_ << ") in the event - Please verify the name of the muon collection";
97  return result;
98  }
99 
100  LogDebug("CosmicRegionalSeedGenerator") << "Muons collection size = " << muonsHandle->size();
101 
102  //get the jet collection
103  edm::Handle<CaloJetCollection> caloJetsHandle;
104  event.getByLabel(recoCaloJetsCollection_,caloJetsHandle);
105 
106  //get the propagator
107  edm::ESHandle<Propagator> thePropagator;
108  es.get<TrackingComponentsRecord>().get(thePropagatorName_, thePropagator); // thePropagatorName = "AnalyticalPropagator"
109 
110  //get tracker geometry
111  edm::ESHandle<TrackerGeometry> theTrackerGeometry;
112  es.get<TrackerDigiGeometryRecord>().get(theTrackerGeometry);
113  //const TrackerGeometry& theTracker(*theTrackerGeometry);
114  DetId outerid;
115 
116 
117  //definition of the region
118  //+++++++++++++++++++++++++
119 
120  int nmuons = 0;
121  for (reco::MuonCollection::const_iterator staMuon = muonsHandle->begin(); staMuon != muonsHandle->end(); ++staMuon) {
122 
123  //select sta muons
124  if (!staMuon->isStandAloneMuon()) {
125  LogDebug("CosmicRegionalSeedGenerator") << "This muon is not a stand alone muon";
126  continue;
127  }
128 
129  //bit 25 as a coverage -1.4 < eta < 1.4
130  if ( abs( staMuon->standAloneMuon()->eta() ) > 1.5 ) continue;
131 
132  //debug
133  nmuons++;
134  LogDebug("CosmicRegionalSeedGenerator") << "Muon stand alone found in the collection - in muons chambers: \n "
135  << "Position = " << staMuon->standAloneMuon()->outerPosition() << "\n "
136  << "Momentum = " << staMuon->standAloneMuon()->outerMomentum() << "\n "
137  << "Eta = " << staMuon->standAloneMuon()->eta() << "\n "
138  << "Phi = " << staMuon->standAloneMuon()->phi();
139 
140  //initial position, momentum, charge
141 
142  GlobalPoint initialRegionPosition(staMuon->standAloneMuon()->referencePoint().x(), staMuon->standAloneMuon()->referencePoint().y(), staMuon->standAloneMuon()->referencePoint().z());
143  GlobalVector initialRegionMomentum(staMuon->standAloneMuon()->momentum().x(), staMuon->standAloneMuon()->momentum().y(), staMuon->standAloneMuon()->momentum().z());
144  int charge = (int) staMuon->standAloneMuon()->charge();
145 
146  LogDebug("CosmicRegionalSeedGenerator") << "Initial region - Reference point of the sta muon: \n "
147  << "Position = " << initialRegionPosition << "\n "
148  << "Momentum = " << initialRegionMomentum << "\n "
149  << "Eta = " << initialRegionPosition.eta() << "\n "
150  << "Phi = " << initialRegionPosition.phi() << "\n "
151  << "Charge = " << charge;
152 
153  //propagation on the last layers of TOB
154  if ( staMuon->standAloneMuon()->outerPosition().y()>0 ) initialRegionMomentum *=-1;
155  GlobalTrajectoryParameters glb_parameters(initialRegionPosition,
156  initialRegionMomentum,
157  charge,
158  thePropagator->magneticField());
159  FreeTrajectoryState fts(glb_parameters);
160  StateOnTrackerBound onBounds(thePropagator.product());
161  TrajectoryStateOnSurface outer = onBounds(fts);
162 
163  if (!outer.isValid())
164  {
165  //edm::LogError("FailedPropagation") << "Trajectory state on surface not valid" ;
166  LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid" ;
167  continue;
168  }
169 
170 
171  //final position & momentum
172  GlobalPoint regionPosition = outer.globalPosition();
173  GlobalVector regionMom = outer.globalMomentum();
174 
175  LogDebug("CosmicRegionalSeedGenerator") << "Region after propagation: \n "
176  << "Position = " << outer.globalPosition() << "\n "
177  << "Momentum = " << outer.globalMomentum() << "\n "
178  << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
179  << "Eta = " << outer.globalPosition().eta() << "\n "
180  << "Phi = " << outer.globalPosition().phi();
181 
182 
183  //step back
184  double stepBack = 1;
185  GlobalPoint center = regionPosition + stepBack * regionMom.unit();
186  GlobalVector v = stepBack * regionMom.unit();
187  LogDebug("CosmicRegionalSeedGenerator") << "Step back vector = " << v << "\n";
188 
189  //exclude region built in jets
190  if ( doJetsExclusionCheck_ ) {
191  double delta_R_min = 1000.;
192  for ( CaloJetCollection::const_iterator jet = caloJetsHandle->begin (); jet != caloJetsHandle->end(); jet++ ) {
193  if ( jet->pt() < jetsPtMin_ ) continue;
194 
195  double deta = center.eta() - jet->eta();
196  double dphi = fabs( center.phi() - jet->phi() );
197  if ( dphi > TMath::Pi() ) dphi = 2*TMath::Pi() - dphi;
198 
199  double delta_R = sqrt(deta*deta + dphi*dphi);
200  if ( delta_R < delta_R_min ) delta_R_min = delta_R;
201 
202  }//end loop on jets
203 
204  if ( delta_R_min < deltaRExclusionSize_ ) {
205  LogDebug("CosmicRegionalSeedGenerator") << "Region built too close from a jet";
206  continue;
207  }
208  }//end if doJetsExclusionCheck
209 
210 
211  //definition of the region
212  CosmicTrackingRegion *etaphiRegion = new CosmicTrackingRegion((-1)*regionMom,
213  center,
214  ptMin_,
215  rVertex_,
216  zVertex_,
217  deltaEta_,
218  deltaPhi_,
219  regionPSet
220  );
221 
222 
223 
224  //return the result
225  result.push_back(etaphiRegion);
226 
227  LogDebug("CosmicRegionalSeedGenerator") << "Final CosmicTrackingRegion \n "
228  << "Position = "<< center << "\n "
229  << "Direction = "<< etaphiRegion->direction() << "\n "
230  << "Distance from the region on the layer = " << (regionPosition -center).mag() << "\n "
231  << "Eta = " << center.eta() << "\n "
232  << "Phi = " << center.phi();
233 
234 
235  }//end loop on muons
236 
237  }//end if SeedOnStaMuon
238 
239 
240 
241 
242 
243  //________________________________________
244  //
245  //Seeding on cosmic muons (MC && Datas)
246  //________________________________________
247 
248 
249  if(regionBase_=="seedOnCosmicMuon") {
250 
251  LogDebug("CosmicRegionalSeedGenerator") << "Seeding on cosmic muons tracks";
252 
253  //get collections
254  //+++++++++++++++
255 
256  //get the muon collection
257  edm::Handle<reco::TrackCollection> cosmicMuonsHandle;
258  event.getByLabel(recoTrackMuonsCollection_,cosmicMuonsHandle);
259  if (!cosmicMuonsHandle.isValid())
260  {
261  edm::LogError("CollectionNotFound") << "Error::No cosmic muons collection (" << recoTrackMuonsCollection_ << ") in the event - Please verify the name of the muon reco track collection";
262  return result;
263  }
264 
265  LogDebug("CosmicRegionalSeedGenerator") << "Cosmic muons tracks collection size = " << cosmicMuonsHandle->size();
266 
267  //get the jet collection
268  edm::Handle<CaloJetCollection> caloJetsHandle;
269  event.getByLabel(recoCaloJetsCollection_,caloJetsHandle);
270 
271  //get the propagator
272  edm::ESHandle<Propagator> thePropagator;
273  es.get<TrackingComponentsRecord>().get(thePropagatorName_, thePropagator); // thePropagatorName = "AnalyticalPropagator"
274 
275  //get tracker geometry
276  edm::ESHandle<TrackerGeometry> theTrackerGeometry;
277  es.get<TrackerDigiGeometryRecord>().get(theTrackerGeometry);
278  DetId outerid;
279 
280 
281  //definition of the region
282  //+++++++++++++++++++++++++
283 
284  int nmuons = 0;
285  for (reco::TrackCollection::const_iterator cosmicMuon = cosmicMuonsHandle->begin(); cosmicMuon != cosmicMuonsHandle->end(); ++cosmicMuon) {
286 
287  //bit 25 as a coverage -1.4 < eta < 1.4
288  if ( abs( cosmicMuon->eta() ) > 1.5 ) continue;
289 
290  nmuons++;
291 
292  //initial position, momentum, charge
293  GlobalPoint initialRegionPosition(cosmicMuon->referencePoint().x(), cosmicMuon->referencePoint().y(), cosmicMuon->referencePoint().z());
294  GlobalVector initialRegionMomentum(cosmicMuon->momentum().x(), cosmicMuon->momentum().y(), cosmicMuon->momentum().z());
295  int charge = (int) cosmicMuon->charge();
296 
297  LogDebug("CosmicRegionalSeedGenerator") << "Position and momentum of the muon track in the muon chambers: \n "
298  << "x = " << cosmicMuon->outerPosition().x() << "\n "
299  << "y = " << cosmicMuon->outerPosition().y() << "\n "
300  << "y = " << cosmicMuon->pt() << "\n "
301  << "Initial region - Reference point of the cosmic muon track: \n "
302  << "Position = " << initialRegionPosition << "\n "
303  << "Momentum = " << initialRegionMomentum << "\n "
304  << "Eta = " << initialRegionPosition.eta() << "\n "
305  << "Phi = " << initialRegionPosition.phi() << "\n "
306  << "Charge = " << charge;
307 
308  //propagation on the last layers of TOB
309  if ( cosmicMuon->outerPosition().y()>0 && cosmicMuon->momentum().y()<0 ) initialRegionMomentum *=-1;
310  GlobalTrajectoryParameters glb_parameters(initialRegionPosition,
311  initialRegionMomentum,
312  charge,
313  thePropagator->magneticField());
314  FreeTrajectoryState fts(glb_parameters);
315  StateOnTrackerBound onBounds(thePropagator.product());
316  TrajectoryStateOnSurface outer = onBounds(fts);
317 
318  if (!outer.isValid())
319  {
320  //edm::LogError("FailedPropagation") << "Trajectory state on surface not valid" ;
321  LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid" ;
322  continue;
323  }
324 
325 
326  //final position & momentum
327  GlobalPoint regionPosition = outer.globalPosition();
328  GlobalVector regionMom = outer.globalMomentum();
329 
330  LogDebug("CosmicRegionalSeedGenerator") << "Region after propagation: \n "
331  << "Position = " << outer.globalPosition() << "\n "
332  << "Momentum = " << outer.globalMomentum() << "\n "
333  << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
334  << "Eta = " << outer.globalPosition().eta() << "\n "
335  << "Phi = " << outer.globalPosition().phi();
336 
337 
338  //step back
339  double stepBack = 1;
340  GlobalPoint center = regionPosition + stepBack * regionMom.unit();
341  GlobalVector v = stepBack * regionMom.unit();
342  LogDebug("CosmicRegionalSeedGenerator") << "Step back vector = " << v << "\n";
343 
344  //exclude region built in jets
345  if ( doJetsExclusionCheck_ ) {
346  double delta_R_min = 1000.;
347  for ( CaloJetCollection::const_iterator jet = caloJetsHandle->begin (); jet != caloJetsHandle->end(); jet++ ) {
348  if ( jet->pt() < jetsPtMin_ ) continue;
349 
350  double deta = center.eta() - jet->eta();
351  double dphi = fabs( center.phi() - jet->phi() );
352  if ( dphi > TMath::Pi() ) dphi = 2*TMath::Pi() - dphi;
353 
354  double delta_R = sqrt(deta*deta + dphi*dphi);
355  if ( delta_R < delta_R_min ) delta_R_min = delta_R;
356 
357  }//end loop on jets
358 
359  if ( delta_R_min < deltaRExclusionSize_ ) {
360  LogDebug("CosmicRegionalSeedGenerator") << "Region built too close from a jet";
361  continue;
362  }
363  }// end if doJetsExclusionCheck
364 
365  //definition of the region
366  CosmicTrackingRegion *etaphiRegion = new CosmicTrackingRegion((-1)*regionMom,
367  center,
368  ptMin_,
369  rVertex_,
370  zVertex_,
371  deltaEta_,
372  deltaPhi_,
373  regionPSet
374  );
375 
376 
377  //return the result
378  result.push_back(etaphiRegion);
379 
380  LogDebug("CosmicRegionalSeedGenerator") << "Final CosmicTrackingRegion \n "
381  << "Position = "<< center << "\n "
382  << "Direction = "<< etaphiRegion->direction() << "\n "
383  << "Distance from the region on the layer = " << (regionPosition -center).mag() << "\n "
384  << "Eta = " << center.eta() << "\n "
385  << "Phi = " << center.phi();
386 
387  }//end loop on muons
388 
389  }//end if SeedOnCosmicMuon
390 
391 
392  //________________________________________
393  //
394  //Seeding on L2 muons (MC && Datas)
395  //________________________________________
396 
397  if(regionBase_=="seedOnL2Muon") {
398 
399  LogDebug("CosmicRegionalSeedGenerator") << "Seeding on L2 muons";
400 
401  //get collections
402  //+++++++++++++++
403 
404  //get the muon collection
406  event.getByLabel(recoL2MuonsCollection_,L2MuonsHandle);
407 
408  if (!L2MuonsHandle.isValid())
409  {
410  edm::LogError("CollectionNotFound") << "Error::No L2 muons collection (" << recoL2MuonsCollection_ <<") in the event - Please verify the name of the L2 muon collection";
411  return result;
412  }
413 
414  LogDebug("CosmicRegionalSeedGenerator") << "L2 muons collection size = " << L2MuonsHandle->size();
415 
416  //get the propagator
417  edm::ESHandle<Propagator> thePropagator;
418  es.get<TrackingComponentsRecord>().get(thePropagatorName_, thePropagator); // thePropagatorName = "AnalyticalPropagator"
419 
420  //get tracker geometry
421  edm::ESHandle<TrackerGeometry> theTrackerGeometry;
422  es.get<TrackerDigiGeometryRecord>().get(theTrackerGeometry);
423  DetId outerid;
424 
425 
426  //definition of the region
427  //+++++++++++++++++++++++++
428 
429  int nmuons = 0;
430  for (reco::RecoChargedCandidateCollection::const_iterator L2Muon = L2MuonsHandle->begin(); L2Muon != L2MuonsHandle->end(); ++L2Muon) {
431  reco::TrackRef tkL2Muon = L2Muon->get<reco::TrackRef>();
432 
433  //bit 25 as a coverage -1.4 < eta < 1.4
434  if ( abs( tkL2Muon->eta() ) > 1.5 ) continue;
435 
436  nmuons++;
437 
438  //initial position, momentum, charge
439  GlobalPoint initialRegionPosition(tkL2Muon->referencePoint().x(), tkL2Muon->referencePoint().y(), tkL2Muon->referencePoint().z());
440  GlobalVector initialRegionMomentum(tkL2Muon->momentum().x(), tkL2Muon->momentum().y(), tkL2Muon->momentum().z());
441  int charge = (int) tkL2Muon->charge();
442 
443  LogDebug("CosmicRegionalSeedGenerator") << "Position and momentum of the L2 muon track in the muon chambers: \n "
444  << "x = " << tkL2Muon->outerPosition().x() << "\n "
445  << "y = " << tkL2Muon->outerPosition().y() << "\n "
446  << "y = " << tkL2Muon->pt() << "\n "
447  << "Initial region - Reference point of the L2 muon track: \n "
448  << "Position = " << initialRegionPosition << "\n "
449  << "Momentum = " << initialRegionMomentum << "\n "
450  << "Eta = " << initialRegionPosition.eta() << "\n "
451  << "Phi = " << initialRegionPosition.phi() << "\n "
452  << "Charge = " << charge;
453 
454 
455  //seeding only in the bottom
456  if ( tkL2Muon->outerPosition().y() > 0 )
457  {
458  LogDebug("CosmicRegionalSeedGenerator") << "L2 muon in the TOP --- Region not created";
459  return result;
460  }
461 
462  GlobalTrajectoryParameters glb_parameters(initialRegionPosition,
463  initialRegionMomentum,
464  charge,
465  thePropagator->magneticField());
466  FreeTrajectoryState fts(glb_parameters);
467  StateOnTrackerBound onBounds(thePropagator.product());
468  TrajectoryStateOnSurface outer = onBounds(fts);
469 
470  if (!outer.isValid())
471  {
472  //edm::LogError("FailedPropagation") << "Trajectory state on surface not valid" ;
473  LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid" ;
474  continue;
475  }
476 
477 
478  //final position & momentum
479  GlobalPoint regionPosition = outer.globalPosition();
480  GlobalVector regionMom = outer.globalMomentum();
481 
482  LogDebug("CosmicRegionalSeedGenerator") << "Region after propagation: \n "
483  << "Position = " << outer.globalPosition() << "\n "
484  << "Momentum = " << outer.globalMomentum() << "\n "
485  << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
486  << "Eta = " << outer.globalPosition().eta() << "\n "
487  << "Phi = " << outer.globalPosition().phi();
488 
489 
490  //step back
491  double stepBack = 1;
492  GlobalPoint center = regionPosition + stepBack * regionMom.unit();
493  GlobalVector v = stepBack * regionMom.unit();
494  LogDebug("CosmicRegionalSeedGenerator") << "Step back vector = " << v << "\n";
495 
496 
497  //definition of the region
498  CosmicTrackingRegion *etaphiRegion = new CosmicTrackingRegion((-1)*regionMom,
499  center,
500  ptMin_,
501  rVertex_,
502  zVertex_,
503  deltaEta_,
504  deltaPhi_,
505  regionPSet
506  );
507 
508  result.push_back(etaphiRegion);
509 
510  LogDebug("CosmicRegionalSeedGenerator") << "Final L2TrackingRegion \n "
511  << "Position = "<< center << "\n "
512  << "Direction = "<< etaphiRegion->direction() << "\n "
513  << "Distance from the region on the layer = " << (regionPosition -center).mag() << "\n "
514  << "Eta = " << center.eta() << "\n "
515  << "Phi = " << center.phi();
516 
517 
518  }//end loop on muons
519 
520  }//end if SeedOnL2Muon
521 
522  return result;
523 
524 }
525 
#define LogDebug(id)
const double Pi
T getParameter(std::string const &) const
T perp() const
Definition: PV3DBase.h:66
Geom::Phi< T > phi() const
Definition: PV3DBase.h:63
#define abs(x)
Definition: mlp_lapack.h:159
virtual GlobalVector direction() const
the direction around which region is constructed
double charge(const std::vector< uint8_t > &Ampls)
CosmicRegionalSeedGenerator(const edm::ParameterSet &conf)
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
T sqrt(T t)
Definition: SSEVec.h:28
T z() const
Definition: PV3DBase.h:58
tuple result
Definition: query.py:137
double delta_R(double eta1, double phi1, double eta2, double phi2)
Definition: AnglesUtil.h:116
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
bool isValid() const
Definition: HandleBase.h:76
tuple conf
Definition: dbtoconf.py:185
Vector3DBase unit() const
Definition: Vector3DBase.h:57
Definition: DetId.h:20
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:62
T eta() const
Definition: PV3DBase.h:70
virtual std::vector< TrackingRegion * > regions(const edm::Event &event, const edm::EventSetup &es) const
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:239
mathSSE::Vec4< T > v