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_=="") {
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  DetId outerid;
109 
110 
111  //definition of the region
112  //+++++++++++++++++++++++++
113 
114  int nmuons = 0;
115  for (reco::MuonCollection::const_iterator staMuon = muonsHandle->begin(); staMuon != muonsHandle->end(); ++staMuon) {
116 
117  //select sta muons
118  if (!staMuon->isStandAloneMuon()) {
119  LogDebug("CosmicRegionalSeedGenerator") << "This muon is not a stand alone muon";
120  continue;
121  }
122 
123  //bit 25 as a coverage -1.4 < eta < 1.4
124  if ( abs( staMuon->standAloneMuon()->eta() ) > 1.5 ) continue;
125 
126  //debug
127  nmuons++;
128  LogDebug("CosmicRegionalSeedGenerator") << "Muon stand alone found in the collection - in muons chambers: \n "
129  << "Position = " << staMuon->standAloneMuon()->outerPosition() << "\n "
130  << "Momentum = " << staMuon->standAloneMuon()->outerMomentum() << "\n "
131  << "Eta = " << staMuon->standAloneMuon()->eta() << "\n "
132  << "Phi = " << staMuon->standAloneMuon()->phi();
133 
134  //initial position, momentum, charge
135 
136  GlobalPoint initialRegionPosition(staMuon->standAloneMuon()->referencePoint().x(), staMuon->standAloneMuon()->referencePoint().y(), staMuon->standAloneMuon()->referencePoint().z());
137  GlobalVector initialRegionMomentum(staMuon->standAloneMuon()->momentum().x(), staMuon->standAloneMuon()->momentum().y(), staMuon->standAloneMuon()->momentum().z());
138  int charge = (int) staMuon->standAloneMuon()->charge();
139 
140  LogDebug("CosmicRegionalSeedGenerator") << "Initial region - Reference point of the sta muon: \n "
141  << "Position = " << initialRegionPosition << "\n "
142  << "Momentum = " << initialRegionMomentum << "\n "
143  << "Eta = " << initialRegionPosition.eta() << "\n "
144  << "Phi = " << initialRegionPosition.phi() << "\n "
145  << "Charge = " << charge;
146 
147  //propagation on the last layers of TOB
148  if ( staMuon->standAloneMuon()->outerPosition().y()>0 ) initialRegionMomentum *=-1;
149  GlobalTrajectoryParameters glb_parameters(initialRegionPosition,
150  initialRegionMomentum,
151  charge,
152  thePropagator->magneticField());
153  FreeTrajectoryState fts(glb_parameters);
154  StateOnTrackerBound onBounds(thePropagator.product());
155  TrajectoryStateOnSurface outer = onBounds(fts);
156 
157  if (!outer.isValid())
158  {
159  //edm::LogError("FailedPropagation") << "Trajectory state on surface not valid" ;
160  LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid" ;
161  continue;
162  }
163 
164 
165  //final position & momentum
166  GlobalPoint regionPosition = outer.globalPosition();
167  GlobalVector regionMom = outer.globalMomentum();
168 
169  LogDebug("CosmicRegionalSeedGenerator") << "Region after propagation: \n "
170  << "Position = " << outer.globalPosition() << "\n "
171  << "Momentum = " << outer.globalMomentum() << "\n "
172  << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
173  << "Eta = " << outer.globalPosition().eta() << "\n "
174  << "Phi = " << outer.globalPosition().phi();
175 
176 
177  //step back
178  double stepBack = 1;
179  GlobalPoint center = regionPosition + stepBack * regionMom.unit();
180  GlobalVector v = stepBack * regionMom.unit();
181  LogDebug("CosmicRegionalSeedGenerator") << "Step back vector = " << v << "\n";
182 
183  //exclude region built in jets
184  if ( doJetsExclusionCheck_ ) {
185  double delta_R_min = 1000.;
186  for ( CaloJetCollection::const_iterator jet = caloJetsHandle->begin (); jet != caloJetsHandle->end(); jet++ ) {
187  if ( jet->pt() < jetsPtMin_ ) continue;
188 
189  double deta = center.eta() - jet->eta();
190  double dphi = fabs( center.phi() - jet->phi() );
191  if ( dphi > TMath::Pi() ) dphi = 2*TMath::Pi() - dphi;
192 
193  double delta_R = sqrt(deta*deta + dphi*dphi);
194  if ( delta_R < delta_R_min ) delta_R_min = delta_R;
195 
196  }//end loop on jets
197 
198  if ( delta_R_min < deltaRExclusionSize_ ) {
199  LogDebug("CosmicRegionalSeedGenerator") << "Region built too close from a jet";
200  continue;
201  }
202  }//end if doJetsExclusionCheck
203 
204 
205  //definition of the region
206 
207  result.push_back(std::make_unique<CosmicTrackingRegion>((-1)*regionMom,
208  center,
209  ptMin_,
210  rVertex_,
211  zVertex_,
212  deltaEta_,
213  deltaPhi_,
214  regionPSet,
215  measurementTracker));
216 
217  LogDebug("CosmicRegionalSeedGenerator") << "Final CosmicTrackingRegion \n "
218  << "Position = "<< center << "\n "
219  << "Direction = "<< result.back()->direction() << "\n "
220  << "Distance from the region on the layer = " << (regionPosition -center).mag() << "\n "
221  << "Eta = " << center.eta() << "\n "
222  << "Phi = " << center.phi();
223 
224 
225  }//end loop on muons
226 
227  }//end if SeedOnStaMuon
228 
229 
230 
231 
232 
233  //________________________________________
234  //
235  //Seeding on cosmic muons (MC && Datas)
236  //________________________________________
237 
238 
239  if(regionBase_=="seedOnCosmicMuon") {
240 
241  LogDebug("CosmicRegionalSeedGenerator") << "Seeding on cosmic muons tracks";
242 
243  //get collections
244  //+++++++++++++++
245 
246  //get the muon collection
247  edm::Handle<reco::TrackCollection> cosmicMuonsHandle;
248  event.getByToken(recoTrackMuonsToken_,cosmicMuonsHandle);
249  if (!cosmicMuonsHandle.isValid())
250  {
251  edm::LogError("CollectionNotFound") << "Error::No cosmic muons collection (" << recoTrackMuonsCollection_ << ") in the event - Please verify the name of the muon reco track collection";
252  return result;
253  }
254 
255  LogDebug("CosmicRegionalSeedGenerator") << "Cosmic muons tracks collection size = " << cosmicMuonsHandle->size();
256 
257  //get the jet collection
258  edm::Handle<CaloJetCollection> caloJetsHandle;
259  event.getByToken(recoCaloJetsToken_,caloJetsHandle);
260 
261  //get the propagator
262  edm::ESHandle<Propagator> thePropagator;
263  es.get<TrackingComponentsRecord>().get(thePropagatorName_, thePropagator); // thePropagatorName = "AnalyticalPropagator"
264 
265  //get tracker geometry
266  edm::ESHandle<TrackerGeometry> theTrackerGeometry;
267  es.get<TrackerDigiGeometryRecord>().get(theTrackerGeometry);
268  DetId outerid;
269 
270 
271  //definition of the region
272  //+++++++++++++++++++++++++
273 
274  int nmuons = 0;
275  for (reco::TrackCollection::const_iterator cosmicMuon = cosmicMuonsHandle->begin(); cosmicMuon != cosmicMuonsHandle->end(); ++cosmicMuon) {
276 
277  //bit 25 as a coverage -1.4 < eta < 1.4
278  if ( abs( cosmicMuon->eta() ) > 1.5 ) continue;
279 
280  nmuons++;
281 
282  //initial position, momentum, charge
283  GlobalPoint initialRegionPosition(cosmicMuon->referencePoint().x(), cosmicMuon->referencePoint().y(), cosmicMuon->referencePoint().z());
284  GlobalVector initialRegionMomentum(cosmicMuon->momentum().x(), cosmicMuon->momentum().y(), cosmicMuon->momentum().z());
285  int charge = (int) cosmicMuon->charge();
286 
287  LogDebug("CosmicRegionalSeedGenerator") << "Position and momentum of the muon track in the muon chambers: \n "
288  << "x = " << cosmicMuon->outerPosition().x() << "\n "
289  << "y = " << cosmicMuon->outerPosition().y() << "\n "
290  << "y = " << cosmicMuon->pt() << "\n "
291  << "Initial region - Reference point of the cosmic muon track: \n "
292  << "Position = " << initialRegionPosition << "\n "
293  << "Momentum = " << initialRegionMomentum << "\n "
294  << "Eta = " << initialRegionPosition.eta() << "\n "
295  << "Phi = " << initialRegionPosition.phi() << "\n "
296  << "Charge = " << charge;
297 
298  //propagation on the last layers of TOB
299  if ( cosmicMuon->outerPosition().y()>0 && cosmicMuon->momentum().y()<0 ) initialRegionMomentum *=-1;
300  GlobalTrajectoryParameters glb_parameters(initialRegionPosition,
301  initialRegionMomentum,
302  charge,
303  thePropagator->magneticField());
304  FreeTrajectoryState fts(glb_parameters);
305  StateOnTrackerBound onBounds(thePropagator.product());
306  TrajectoryStateOnSurface outer = onBounds(fts);
307 
308  if (!outer.isValid())
309  {
310  //edm::LogError("FailedPropagation") << "Trajectory state on surface not valid" ;
311  LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid" ;
312  continue;
313  }
314 
315 
316  //final position & momentum
317  GlobalPoint regionPosition = outer.globalPosition();
318  GlobalVector regionMom = outer.globalMomentum();
319 
320  LogDebug("CosmicRegionalSeedGenerator") << "Region after propagation: \n "
321  << "Position = " << outer.globalPosition() << "\n "
322  << "Momentum = " << outer.globalMomentum() << "\n "
323  << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
324  << "Eta = " << outer.globalPosition().eta() << "\n "
325  << "Phi = " << outer.globalPosition().phi();
326 
327 
328  //step back
329  double stepBack = 1;
330  GlobalPoint center = regionPosition + stepBack * regionMom.unit();
331  GlobalVector v = stepBack * regionMom.unit();
332  LogDebug("CosmicRegionalSeedGenerator") << "Step back vector = " << v << "\n";
333 
334  //exclude region built in jets
335  if ( doJetsExclusionCheck_ ) {
336  double delta_R_min = 1000.;
337  for ( CaloJetCollection::const_iterator jet = caloJetsHandle->begin (); jet != caloJetsHandle->end(); jet++ ) {
338  if ( jet->pt() < jetsPtMin_ ) continue;
339 
340  double deta = center.eta() - jet->eta();
341  double dphi = fabs( center.phi() - jet->phi() );
342  if ( dphi > TMath::Pi() ) dphi = 2*TMath::Pi() - dphi;
343 
344  double delta_R = sqrt(deta*deta + dphi*dphi);
345  if ( delta_R < delta_R_min ) delta_R_min = delta_R;
346 
347  }//end loop on jets
348 
349  if ( delta_R_min < deltaRExclusionSize_ ) {
350  LogDebug("CosmicRegionalSeedGenerator") << "Region built too close from a jet";
351  continue;
352  }
353  }// end if doJetsExclusionCheck
354 
355  //definition of the region
356  result.push_back(std::make_unique<CosmicTrackingRegion>((-1)*regionMom,
357  center,
358  ptMin_,
359  rVertex_,
360  zVertex_,
361  deltaEta_,
362  deltaPhi_,
363  regionPSet,
364  measurementTracker));
365 
366  LogDebug("CosmicRegionalSeedGenerator") << "Final CosmicTrackingRegion \n "
367  << "Position = "<< center << "\n "
368  << "Direction = "<< result.back()->direction() << "\n "
369  << "Distance from the region on the layer = " << (regionPosition -center).mag() << "\n "
370  << "Eta = " << center.eta() << "\n "
371  << "Phi = " << center.phi();
372 
373  }//end loop on muons
374 
375  }//end if SeedOnCosmicMuon
376 
377 
378  //________________________________________
379  //
380  //Seeding on L2 muons (MC && Datas)
381  //________________________________________
382 
383  if(regionBase_=="seedOnL2Muon") {
384 
385  LogDebug("CosmicRegionalSeedGenerator") << "Seeding on L2 muons";
386 
387  //get collections
388  //+++++++++++++++
389 
390  //get the muon collection
392  event.getByToken(recoL2MuonsToken_,L2MuonsHandle);
393 
394  if (!L2MuonsHandle.isValid())
395  {
396  edm::LogError("CollectionNotFound") << "Error::No L2 muons collection (" << recoL2MuonsCollection_ <<") in the event - Please verify the name of the L2 muon collection";
397  return result;
398  }
399 
400  LogDebug("CosmicRegionalSeedGenerator") << "L2 muons collection size = " << L2MuonsHandle->size();
401 
402  //get the propagator
403  edm::ESHandle<Propagator> thePropagator;
404  es.get<TrackingComponentsRecord>().get(thePropagatorName_, thePropagator); // thePropagatorName = "AnalyticalPropagator"
405 
406  //get tracker geometry
407  edm::ESHandle<TrackerGeometry> theTrackerGeometry;
408  es.get<TrackerDigiGeometryRecord>().get(theTrackerGeometry);
409  DetId outerid;
410 
411 
412  //definition of the region
413  //+++++++++++++++++++++++++
414 
415  int nmuons = 0;
416  for (reco::RecoChargedCandidateCollection::const_iterator L2Muon = L2MuonsHandle->begin(); L2Muon != L2MuonsHandle->end(); ++L2Muon) {
417  reco::TrackRef tkL2Muon = L2Muon->get<reco::TrackRef>();
418 
419  //bit 25 as a coverage -1.4 < eta < 1.4
420  if ( abs( tkL2Muon->eta() ) > 1.5 ) continue;
421 
422  nmuons++;
423 
424  //initial position, momentum, charge
425  GlobalPoint initialRegionPosition(tkL2Muon->referencePoint().x(), tkL2Muon->referencePoint().y(), tkL2Muon->referencePoint().z());
426  GlobalVector initialRegionMomentum(tkL2Muon->momentum().x(), tkL2Muon->momentum().y(), tkL2Muon->momentum().z());
427  int charge = (int) tkL2Muon->charge();
428 
429  LogDebug("CosmicRegionalSeedGenerator") << "Position and momentum of the L2 muon track in the muon chambers: \n "
430  << "x = " << tkL2Muon->outerPosition().x() << "\n "
431  << "y = " << tkL2Muon->outerPosition().y() << "\n "
432  << "y = " << tkL2Muon->pt() << "\n "
433  << "Initial region - Reference point of the L2 muon track: \n "
434  << "Position = " << initialRegionPosition << "\n "
435  << "Momentum = " << initialRegionMomentum << "\n "
436  << "Eta = " << initialRegionPosition.eta() << "\n "
437  << "Phi = " << initialRegionPosition.phi() << "\n "
438  << "Charge = " << charge;
439 
440 
441  //seeding only in the bottom
442  if ( tkL2Muon->outerPosition().y() > 0 )
443  {
444  LogDebug("CosmicRegionalSeedGenerator") << "L2 muon in the TOP --- Region not created";
445  return result;
446  }
447 
448  GlobalTrajectoryParameters glb_parameters(initialRegionPosition,
449  initialRegionMomentum,
450  charge,
451  thePropagator->magneticField());
452  FreeTrajectoryState fts(glb_parameters);
453  StateOnTrackerBound onBounds(thePropagator.product());
454  TrajectoryStateOnSurface outer = onBounds(fts);
455 
456  if (!outer.isValid())
457  {
458  //edm::LogError("FailedPropagation") << "Trajectory state on surface not valid" ;
459  LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid" ;
460  continue;
461  }
462 
463 
464  //final position & momentum
465  GlobalPoint regionPosition = outer.globalPosition();
466  GlobalVector regionMom = outer.globalMomentum();
467 
468  LogDebug("CosmicRegionalSeedGenerator") << "Region after propagation: \n "
469  << "Position = " << outer.globalPosition() << "\n "
470  << "Momentum = " << outer.globalMomentum() << "\n "
471  << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
472  << "Eta = " << outer.globalPosition().eta() << "\n "
473  << "Phi = " << outer.globalPosition().phi();
474 
475 
476  //step back
477  double stepBack = 1;
478  GlobalPoint center = regionPosition + stepBack * regionMom.unit();
479  GlobalVector v = stepBack * regionMom.unit();
480  LogDebug("CosmicRegionalSeedGenerator") << "Step back vector = " << v << "\n";
481 
482 
483  //definition of the region
484  result.push_back(std::make_unique<CosmicTrackingRegion>((-1)*regionMom,
485  center,
486  ptMin_,
487  rVertex_,
488  zVertex_,
489  deltaEta_,
490  deltaPhi_,
491  regionPSet,
492  measurementTracker));
493 
494  LogDebug("CosmicRegionalSeedGenerator") << "Final L2TrackingRegion \n "
495  << "Position = "<< center << "\n "
496  << "Direction = "<< result.back()->direction() << "\n "
497  << "Distance from the region on the layer = " << (regionPosition -center).mag() << "\n "
498  << "Eta = " << center.eta() << "\n "
499  << "Phi = " << center.phi();
500 
501 
502  }//end loop on muons
503 
504  }//end if SeedOnL2Muon
505 
506  return result;
507 
508 }
509 
#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:14
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:245
edm::EDGetTokenT< reco::MuonCollection > recoMuonsToken_
double delta_R(double eta1, double phi1, double eta2, double phi2)
Definition: AnglesUtil.h:116
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
Definition: DetId.h:18
T const * product() const
Definition: Handle.h:81
const T & get() const
Definition: EventSetup.h:55
T eta() const
Definition: PV3DBase.h:76
fixed size matrix
HLT enums.
bool isUninitialized() const
Definition: EDGetToken.h:73
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