CMS 3D CMS Logo

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