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