CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonShowerInformationFiller.cc
Go to the documentation of this file.
1 
13 
14 // system include files
15 #include <memory>
16 #include <algorithm>
17 #include <iostream>
18 
19 // user include files
52 
55 
61 
62 using namespace std;
63 using namespace edm;
64 
65 //
66 // Constructor
67 //
69  theService(0),
70  theDTRecHitLabel(par.getParameter<InputTag>("DTRecSegmentLabel")),
71  theCSCRecHitLabel(par.getParameter<InputTag>("CSCRecSegmentLabel")),
72  theCSCSegmentsLabel(par.getParameter<InputTag>("CSCSegmentLabel")),
73  theDT4DRecSegmentLabel(par.getParameter<InputTag>("DT4DRecSegmentLabel"))
74 {
75 
80 
81 
82  edm::ParameterSet serviceParameters = par.getParameter<edm::ParameterSet>("ServiceParameters");
83  theService = new MuonServiceProxy(serviceParameters);
84 
85  theTrackerRecHitBuilderName = par.getParameter<string>("TrackerRecHitBuilder");
86  theMuonRecHitBuilderName = par.getParameter<string>("MuonRecHitBuilder");
87 
88  theCacheId_TRH = 0;
89  theCacheId_MT = 0;
90 
91  category_ = "MuonShowerInformationFiller";
92 
93  for (int istat = 0; istat < 4; istat++) {
94  theStationShowerDeltaR.push_back(0.);
95  theStationShowerTSize.push_back(0.);
96  theAllStationHits.push_back(0);
97  theCorrelatedStationHits.push_back(0);
98  }
99 
100 }
101 
102 //
103 // Destructor
104 //
106  if (theService) delete theService;
107 }
108 
109 //
110 //Fill the MuonShower struct
111 //
113 
114  reco::MuonShower returnShower;
115 
116  // Update the services
117  theService->update(iSetup);
118  setEvent(iEvent);
119  setServices(theService->eventSetup());
120 
121  fillHitsByStation(muon);
122  std::vector<int> nStationHits = theAllStationHits;
123  std::vector<int> nStationCorrelatedHits = theCorrelatedStationHits;
124  std::vector<float> stationShowerSizeT = theStationShowerTSize;
125  std::vector<float> stationShowerDeltaR = theStationShowerDeltaR;
126 
127  returnShower.nStationHits = nStationHits;
128  returnShower.nStationCorrelatedHits = nStationCorrelatedHits;
129  returnShower.stationShowerSizeT = stationShowerSizeT;
130  returnShower.stationShowerDeltaR = stationShowerDeltaR;
131 
132  return returnShower;
133 
134 }
135 
136 //
137 // Set Event
138 //
140 
141  // get all the necesary products
142  event.getByToken(theDTRecHitToken, theDTRecHits);
143  event.getByToken(theCSCRecHitToken, theCSCRecHits);
144  event.getByToken(theCSCSegmentsToken, theCSCSegments);
145  event.getByToken(theDT4DRecSegmentToken, theDT4DRecSegments);
146 
147  for (int istat = 0; istat < 4; istat++) {
148  theStationShowerDeltaR.at(istat) = 0.;
149  theStationShowerTSize.at(istat) = 0.;
150  theAllStationHits.at(istat) = 0;
151  theCorrelatedStationHits.at(istat) = 0;
152  }
153 }
154 
155 
156 //
157 // Set services
158 //
160 
161  // DetLayer Geometry
163  setup.get<IdealMagneticFieldRecord>().get(theField);
165  setup.get<MuonGeometryRecord>().get(theCSCGeometry);
166  setup.get<MuonGeometryRecord>().get(theDTGeometry);
167 
168  // Transient Rechit Builders
169  unsigned long long newCacheId_TRH = setup.get<TransientRecHitRecord>().cacheIdentifier();
170  if ( newCacheId_TRH != theCacheId_TRH ) {
173  }
174 
175 }
176 
177 //
178 // Get hits owned by segments
179 //
184 
186 
187  DetId geoId = geomDet->geographicalId();
188 
189  if (geoId.subdetId() == MuonSubdetId::DT) {
190 
191  DTChamberId chamberId(geoId.rawId());
192 
193  // loop on segments 4D
194  DTRecSegment4DCollection::id_iterator chamberIdIt;
195  for (chamberIdIt = dtSegments->id_begin();
196  chamberIdIt != dtSegments->id_end();
197  ++chamberIdIt){
198 
199  if (*chamberIdIt != chamberId) continue;
200 
201  // Get the range for the corresponding ChamberId
202  DTRecSegment4DCollection::range range = dtSegments->get((*chamberIdIt));
203 
204  for (DTRecSegment4DCollection::const_iterator iseg = range.first;
205  iseg!=range.second;++iseg) {
206  if (iseg->dimension() != 4) continue;
207  segments.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet,&*iseg));
208  }
209  }
210  }
211  else if (geoId.subdetId() == MuonSubdetId::CSC) {
212 
213  CSCDetId did(geoId.rawId());
214 
215  for ( CSCSegmentCollection::id_iterator chamberId = cscSegments->id_begin();
216  chamberId != cscSegments->id_end(); ++chamberId) {
217 
218  if ((*chamberId).chamber() != did.chamber()) continue;
219 
220  // Get the range for the corresponding ChamberId
221  CSCSegmentCollection::range range = cscSegments->get((*chamberId));
222 
223  for (CSCSegmentCollection::const_iterator iseg = range.first;
224  iseg!=range.second;++iseg) {
225  if (iseg->dimension() != 3) continue;
226  segments.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet,&*iseg));
227  }
228  }
229  }
230  else {
231  LogTrace(category_) << "Segments are not built in RPCs" << endl;
232  }
233 
235 
236  if (segments.empty()) return allhitscorrelated;
237 
238  TransientTrackingRecHit::ConstRecHitPointer muonRecHit(segments.front());
239  allhitscorrelated = MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit,2);
240 
241  if (segments.size() == 1) return allhitscorrelated;
242 
243  for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator iseg = segments.begin() + 1;
244  iseg != segments.end(); ++iseg) {
245 
248 
249  for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit1 = hits1.begin();
250  ihit1 != hits1.end(); ++ihit1 ) {
251 
252  bool usedbefore = false;
253  //unused DetId thisID = (*ihit1)->geographicalId();
254  //LocalPoint lp1dinsegHit = (*ihit1)->localPosition();
255  GlobalPoint gp1dinsegHit = (*ihit1)->globalPosition();
256 
257  for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit2 = allhitscorrelated.begin();
258  ihit2 != allhitscorrelated.end(); ++ihit2 ) {
259 
260  //unused DetId thisID2 = (*ihit2)->geographicalId();
261  //LocalPoint lp1dinsegHit2 = (*ihit2)->localPosition();
262  GlobalPoint gp1dinsegHit2 = (*ihit2)->globalPosition();
263 
264  if ( (gp1dinsegHit2 - gp1dinsegHit).mag() < 1.0 ) usedbefore = true;
265 
266  }
267  if ( !usedbefore ) allhitscorrelated.push_back(*ihit1);
268  }
269  }
270 
271  return allhitscorrelated;
272 
273 }
274 
275 
276 //
277 // Find cluster
278 //
281  const GlobalPoint& refpoint) const {
282 
283  if ( muonRecHits.empty() ) return muonRecHits;
284 
285  //clustering step by phi
286  float step = 0.05;
288 
289  stable_sort(muonRecHits.begin(), muonRecHits.end(), AbsLessDPhi(refpoint));
290 
291  for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ihit = muonRecHits.begin(); ihit != muonRecHits.end() - 1; ++ihit) {
292  if (fabs(deltaPhi((*(ihit+1))->globalPosition().phi(), (*ihit)->globalPosition().phi() )) < step) {
293  result.push_back(*ihit);
294  } else {
295  break;
296  }
297  }
298 
299  LogTrace(category_) << "phi front: " << muonRecHits.front()->globalPosition().phi() << endl;
300  LogTrace(category_) << "phi back: " << muonRecHits.back()->globalPosition().phi() << endl;
301 
302  return result;
303 
304 }
305 
306 //
307 //
308 //
311  const GlobalPoint& refpoint) const {
312 
313  if ( muonRecHits.empty() ) return muonRecHits;
314 
315  //clustering step by theta
316  float step = 0.05;
318 
319  stable_sort(muonRecHits.begin(), muonRecHits.end(), AbsLessDTheta(refpoint));
320 
321  for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit = muonRecHits.begin(); ihit != muonRecHits.end() - 1; ++ihit) {
322  if (fabs((*(ihit+1))->globalPosition().theta() - (*ihit)->globalPosition().theta() ) < step) {
323  result.push_back(*ihit);
324  } else {
325  break;
326  }
327  }
328 
329  return result;
330 
331 }
332 
333 //
334 //Used to treat overlap region
335 //
338 
339  if ( muonRecHits.empty() ) return muonRecHits;
340 
341  stable_sort(muonRecHits.begin(), muonRecHits.end(), LessPerp());
342 
343  MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator
344  seedhit = min_element(muonRecHits.begin(), muonRecHits.end(), LessPerp());
345 
346  MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ihigh = seedhit;
347  MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ilow = seedhit;
348 
349  float step = 0.1;
350  while (ihigh != muonRecHits.end()-1 && ( fabs((*(ihigh+1))->globalPosition().perp() - (*ihigh)->globalPosition().perp() ) < step) ) {
351  ihigh++;
352  }
353  while (ilow != muonRecHits.begin() && ( fabs((*ilow)->globalPosition().perp() - (*(ilow -1))->globalPosition().perp()) < step ) ) {
354  ilow--;
355  }
356 
358 
359  return result;
360 
361 }
362 
363 //
364 // Get compatible dets
365 //
366 vector<const GeomDet*> MuonShowerInformationFiller::getCompatibleDets(const reco::Track& track) const {
367 
368  vector<const GeomDet*> total;
369  total.reserve(1000);
370 
371  LogTrace(category_) << "Consider a track " << track.p() << " eta: " << track.eta() << " phi " << track.phi() << endl;
372 
373 
374  TrajectoryStateOnSurface innerTsos = trajectoryStateTransform::innerStateOnSurface(track, *theService->trackingGeometry(), &*theService->magneticField());
375  TrajectoryStateOnSurface outerTsos = trajectoryStateTransform::outerStateOnSurface(track, *theService->trackingGeometry(), &*theService->magneticField());
376 
377  GlobalPoint innerPos = innerTsos.globalPosition();
378  GlobalPoint outerPos = outerTsos.globalPosition();
379 
380  vector<GlobalPoint> allCrossingPoints;
381 
382  const vector<const DetLayer*>& dtlayers = theService->detLayerGeometry()->allDTLayers();
383 
384  for (auto iLayer = dtlayers.begin(); iLayer != dtlayers.end(); ++iLayer) {
385 
386  // crossing points of track with cylinder
387  GlobalPoint xPoint = crossingPoint(innerPos, outerPos, dynamic_cast<const BarrelDetLayer*>(*iLayer));
388 
389  // check if point is inside the detector
390  if ((fabs(xPoint.y()) < 1000.0) && (fabs(xPoint.z()) < 1500 ) &&
391  (!(xPoint.y() == 0 && xPoint.x() == 0 && xPoint.z() == 0))) allCrossingPoints.push_back(xPoint);
392  }
393 
394  stable_sort(allCrossingPoints.begin(), allCrossingPoints.end(), LessMag(innerPos) );
395 
396  vector<const GeomDet*> tempDT;
397 
398  for (vector<GlobalPoint>::const_iterator ipos = allCrossingPoints.begin(); ipos != allCrossingPoints.end(); ++ipos) {
399 
400  tempDT = dtPositionToDets(*ipos);
401  vector<const GeomDet*>::const_iterator begin = tempDT.begin();
402  vector<const GeomDet*>::const_iterator end = tempDT.end();
403 
404  for (; begin!=end;++begin) {
405  total.push_back(*begin);
406  }
407 
408  }
409  allCrossingPoints.clear();
410 
411  const vector<const DetLayer*>& csclayers = theService->detLayerGeometry()->allCSCLayers();
412  for (auto iLayer = csclayers.begin(); iLayer != csclayers.end(); ++iLayer) {
413 
414  GlobalPoint xPoint = crossingPoint(innerPos, outerPos, dynamic_cast<const ForwardDetLayer*>(*iLayer));
415 
416  // check if point is inside the detector
417  if ((fabs(xPoint.y()) < 1000.0) && (fabs(xPoint.z()) < 1500.0)
418  && (!(xPoint.y() == 0 && xPoint.x() == 0 && xPoint.z() == 0))) allCrossingPoints.push_back(xPoint);
419  }
420  stable_sort(allCrossingPoints.begin(), allCrossingPoints.end(), LessMag(innerPos) );
421 
422  vector<const GeomDet*> tempCSC;
423  for (vector<GlobalPoint>::const_iterator ipos = allCrossingPoints.begin(); ipos != allCrossingPoints.end(); ++ipos) {
424 
425  tempCSC = cscPositionToDets(*ipos);
426  vector<const GeomDet*>::const_iterator begin = tempCSC.begin();
427  vector<const GeomDet*>::const_iterator end = tempCSC.end();
428 
429  for (; begin!=end;++begin) {
430  total.push_back(*begin);
431  }
432 
433  }
434 
435  return total;
436 
437 }
438 
439 
440 //
441 // Intersection point of track with barrel layer
442 //
444  const GlobalPoint& p2,
445  const BarrelDetLayer* dl) const {
446 
447  const BoundCylinder& bc = dl->specificSurface();
448  return crossingPoint(p1, p2, bc);
449 
450 }
451 
453  const GlobalPoint& p2,
454  const Cylinder& cyl) const {
455 
456  float radius = cyl.radius();
457 
458  GlobalVector dp = p1 - p2;
459  float slope = dp.x()/dp.y();
460  float a = p1.x() - slope * p1.y();
461 
462  float n2 = (1 + slope * slope);
463  float n1 = 2*a*slope;
464  float n0 = a*a - radius*radius;
465 
466  float y1 = 9999;
467  float y2 = 9999;
468  if ( n1*n1 - 4*n2*n0 > 0 ) {
469  y1 = (-n1 + sqrt(n1*n1 - 4*n2*n0) ) / (2 * n2);
470  y2 = (-n1 - sqrt(n1*n1 - 4*n2*n0) ) / (2 * n2);
471  }
472 
473  float x1 = p1.x() + slope * (y1 - p1.y());
474  float x2 = p1.x() + slope * (y2 - p1.y());
475 
476  float slopeZ = dp.z()/dp.y();
477 
478  float z1 = p1.z() + slopeZ * (y1 - p1.y());
479  float z2 = p1.z() + slopeZ * (y2 - p1.y());
480 
481  // there are two crossing points, return the one that is in the same quadrant as point of extrapolation
482  if ((p2.x()*x1 > 0) && (y1*p2.y() > 0) && (z1*p2.z() > 0)) {
483  return GlobalPoint(x1, y1, z1);
484  } else {
485  return GlobalPoint(x2, y2, z2);
486  }
487 
488 }
489 
490 
491 //
492 // Intersection point of track with a forward layer
493 //
495  const GlobalPoint& p2,
496  const ForwardDetLayer* dl) const {
497 
498  const BoundDisk& bc = dl->specificSurface();
499  return crossingPoint(p1, p2, bc);
500 
501 }
502 
504  const GlobalPoint& p2,
505  const BoundDisk& disk) const {
506 
507  float diskZ = disk.position().z();
508  int endcap = diskZ > 0 ? 1 : (diskZ < 0 ? -1 : 0);
509  diskZ = diskZ + endcap*dynamic_cast<const SimpleDiskBounds&>(disk.bounds()).thickness()/2.;
510 
511  GlobalVector dp = p1 - p2;
512 
513  float slopeZ = dp.z()/dp.y();
514  float y1 = diskZ / slopeZ;
515 
516  float slopeX = dp.z()/dp.x();
517  float x1 = diskZ / slopeX;
518 
519  float z1 = diskZ;
520 
521  if (p2.z()*z1 > 0) {
522  return GlobalPoint(x1, y1, z1);
523  } else {
524  return GlobalPoint(0, 0, 0);
525  }
526 
527 }
528 
529 
530 //
531 // GeomDets along the track in DT
532 //
533 vector<const GeomDet*> MuonShowerInformationFiller::dtPositionToDets(const GlobalPoint& gp) const {
534 
535  int minwheel = -3;
536  int maxwheel = -3;
537  if ( gp.z() < -680.0 ) { minwheel = -3; maxwheel = -3;}
538  else if ( gp.z() < -396.0 ) { minwheel = -2; maxwheel = -1;}
539  else if (gp.z() < -126.8) { minwheel = -2; maxwheel = 0; }
540  else if (gp.z() < 126.8) { minwheel = -1; maxwheel = 1; }
541  else if (gp.z() < 396.0) { minwheel = 0; maxwheel = 2; }
542  else if (gp.z() < 680.0) { minwheel = 1; maxwheel = 2; }
543  else { minwheel = 3; maxwheel = 3; }
544 
545  int station = 5;
546  if ( gp.perp() > 680.0 && gp.perp() < 755.0 ) station = 4;
547  else if ( gp.perp() > 580.0 ) station = 3;
548  else if ( gp.perp() > 480.0 ) station = 2;
549  else if ( gp.perp() > 380.0 ) station = 1;
550  else station = 0;
551 
552  vector<int> sectors;
553 
554  float phistep = M_PI/6;
555 
556  float phigp = (float)gp.phi();
557 
558  if ( fabs(deltaPhi(phigp, 0*phistep)) < phistep ) sectors.push_back(1);
559  if ( fabs(deltaPhi(phigp, phistep)) < phistep ) sectors.push_back(2);
560  if ( fabs(deltaPhi(phigp, 2*phistep)) < phistep ) sectors.push_back(3);
561  if ( fabs(deltaPhi(phigp, 3*phistep)) < phistep ) {
562  sectors.push_back(4);
563  if (station == 4) sectors.push_back(13);
564  }
565  if ( fabs(deltaPhi(phigp, 4*phistep)) < phistep ) sectors.push_back(5);
566  if ( fabs(deltaPhi(phigp, 5*phistep)) < phistep ) sectors.push_back(6);
567  if ( fabs(deltaPhi(phigp, 6*phistep)) < phistep ) sectors.push_back(7);
568  if ( fabs(deltaPhi(phigp, 7*phistep)) < phistep ) sectors.push_back(8);
569  if ( fabs(deltaPhi(phigp, 8*phistep)) < phistep ) sectors.push_back(9);
570  if ( fabs(deltaPhi(phigp, 9*phistep)) < phistep ) {
571  sectors.push_back(10);
572  if (station == 4) sectors.push_back(14);
573  }
574  if ( fabs(deltaPhi(phigp, 10*phistep)) < phistep ) sectors.push_back(11);
575  if ( fabs(deltaPhi(phigp, 11*phistep)) < phistep ) sectors.push_back(12);
576 
577  LogTrace(category_) << "DT position to dets" << endl;
578  LogTrace(category_) << "number of sectors to consider: " << sectors.size() << endl;
579  LogTrace(category_) << "station: " << station << " wheels: " << minwheel << " " << maxwheel << endl;
580 
581  vector<const GeomDet*> result;
582  if (station > 4 || station < 1) return result;
583  if (minwheel > 2 || maxwheel < -2) return result;
584 
585  for (vector<int>::const_iterator isector = sectors.begin(); isector != sectors.end(); ++isector ) {
586  for (int iwheel = minwheel; iwheel != maxwheel + 1; ++iwheel) {
587  DTChamberId chamberid(iwheel, station, (*isector));
588  result.push_back(theService->trackingGeometry()->idToDet(chamberid));
589  }
590  }
591 
592  LogTrace(category_) << "number of GeomDets for this track: " << result.size() << endl;
593 
594  return result;
595 
596 }
597 
598 
599 //
600 // GeomDets along the track in CSC
601 //
602 vector<const GeomDet*> MuonShowerInformationFiller::cscPositionToDets(const GlobalPoint& gp) const {
603 
604  // determine the endcap side
605  int endcap = 0;
606  if (gp.z() > 0) {endcap = 1;} else {endcap = 2;}
607 
608  // determine the csc station and range of rings
609  int station = 5;
610 
611  // check all rings in a station
612  if ( fabs(gp.z()) > 1000. && fabs(gp.z()) < 1055.0 ) {
613  station = 4;
614  }
615  else if ( fabs(gp.z()) > 910.0 && fabs(gp.z()) < 965.0) {
616  station = 3;
617  }
618  else if ( fabs(gp.z()) > 800.0 && fabs(gp.z()) < 860.0) {
619  station = 2;
620  }
621  else if ( fabs(gp.z()) > 570.0 && fabs(gp.z()) < 730.0) {
622  station = 1;
623  }
624 
625  vector<int> sectors;
626 
627  float phistep1 = M_PI/18.; //for all the rings except first rings for stations > 1
628  float phistep2 = M_PI/9.;
629  float phigp = (float)gp.phi();
630 
631  int ring = -1;
632 
633  // determine the ring
634  if (station == 1) {
635 
636 //FIX ME!!!
637  if (gp.perp() > 100 && gp.perp() < 270) ring = 1;
638  else if (gp.perp() > 270 && gp.perp() < 450) ring = 2;
639  else if (gp.perp() > 450 && gp.perp() < 695) ring = 3;
640  else if (gp.perp() > 100 && gp.perp() < 270) ring = 4;
641 
642  }
643  else if (station == 2) {
644 
645  if (gp.perp() > 140 && gp.perp() < 350) ring = 1;
646  else if (gp.perp() > 350 && gp.perp() < 700) ring = 2;
647 
648  }
649  else if (station == 3) {
650 
651  if (gp.perp() > 160 && gp.perp() < 350) ring = 1;
652  else if (gp.perp() > 350 && gp.perp() < 700) ring = 2;
653 
654  }
655  else if (station == 4) {
656 
657  if (gp.perp() > 175 && gp.perp() < 350) ring = 1;
658  else if (gp.perp() > 350 && gp.perp() < 700) ring = 2;
659 
660  }
661 
662  if (station > 1 && ring == 1) {
663 
664  // we have 18 sectors in that case
665  for (int i = 0; i < 18; i++) {
666  if ( fabs(deltaPhi(phigp, i*phistep2)) < phistep2 ) sectors.push_back(i+1);
667  }
668 
669  } else {
670 
671  // we have 36 sectors in that case
672  for (int i = 0; i < 36; i++) {
673  if ( fabs(deltaPhi(phigp, i*phistep1)) < phistep1 ) sectors.push_back(i+1);
674  }
675 }
676 
677  LogTrace(category_) << "CSC position to dets" << endl;
678  LogTrace(category_) << "ring: " << ring << endl;
679  LogTrace(category_) << "endcap: " << endcap << endl;
680  LogTrace(category_) << "station: " << station << endl;
681  LogTrace(category_) << "CSC number of sectors to consider: " << sectors.size() << endl;
682 
683 
684  // check exceptional cases
685  vector<const GeomDet*> result;
686  if (station > 4 || station < 1) return result;
687  if (endcap == 0) return result;
688  if (ring == -1) return result;
689 
690  int minlayer = 1;
691  int maxlayer = 6;
692 
693  for (vector<int>::const_iterator isector = sectors.begin(); isector != sectors.end(); ++isector) {
694  for (int ilayer = minlayer; ilayer != maxlayer + 1; ++ ilayer) {
695  CSCDetId cscid(endcap, station, ring, (*isector), ilayer);
696  result.push_back(theService->trackingGeometry()->idToDet(cscid));
697  }
698  }
699 
700  return result;
701 
702 }
703 
704 //
705 //Set class members
706 //
708 
709  reco::TrackRef track;
710  if ( muon.isGlobalMuon() ) track = muon.globalTrack();
711  else if ( muon.isStandAloneMuon() ) track = muon.outerTrack();
712  else return;
713 
714  // split 1D rechits by station
715  vector<MuonRecHitContainer> muonRecHits(4);
716 
717  // split rechits from segs by station
718  vector<TransientTrackingRecHit::ConstRecHitContainer> muonCorrelatedHits(4);
719 
720  // get vector of GeomDets compatible with a track
721  vector<const GeomDet*> compatibleLayers = getCompatibleDets(*track);
722 
723  // for special cases: CSC station 1
724  MuonRecHitContainer tmpCSC1;
725  bool dtOverlapToCheck = false;
726  bool cscOverlapToCheck = false;
727 
728  for (vector<const GeomDet*>::const_iterator igd = compatibleLayers.begin(); igd != compatibleLayers.end(); igd++ ) {
729 
730  // get det id
731  DetId geoId = (*igd)->geographicalId();
732 
733  // skip tracker hits
734  if (geoId.det()!= DetId::Muon) continue;
735 
736  // DT
737  if ( geoId.subdetId() == MuonSubdetId::DT ) {
738 
739  // get station
740  DTChamberId detid(geoId.rawId());
741  int station = detid.station();
742  int wheel = detid.wheel();
743 
744  // get rechits from segments per station
746  TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_begin = muonCorrelatedHitsTmp.begin();
747  TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_end = muonCorrelatedHitsTmp.end();
748 
749  for (; hits_begin!= hits_end;++hits_begin) {
750  muonCorrelatedHits.at(station-1).push_back(*hits_begin);
751  }
752 
753  //check overlap certain wheels and stations
754  if (abs(wheel) == 2 && station != 4 && station != 1) dtOverlapToCheck = true;
755 
756  // loop over all superlayers of a DT chamber
757  for (int isuperlayer = DTChamberId::minSuperLayerId; isuperlayer != DTChamberId::maxSuperLayerId + 1; ++isuperlayer) {
758  // loop over all layers inside the superlayer
759  for (int ilayer = DTChamberId::minLayerId; ilayer != DTChamberId::maxLayerId+1; ++ilayer) {
760  DTLayerId lid(detid, isuperlayer, ilayer);
761  DTRecHitCollection::range dRecHits = theDTRecHits->get(lid);
762  for (DTRecHitCollection::const_iterator rechit = dRecHits.first; rechit != dRecHits.second;++rechit) {
763  vector<const TrackingRecHit*> subrechits = (*rechit).recHits();
764  for (vector<const TrackingRecHit*>::iterator irechit = subrechits.begin(); irechit != subrechits.end(); ++irechit) {
765  muonRecHits.at(station-1).push_back(MuonTransientTrackingRecHit::specificBuild((&**igd),&**irechit));
766  }
767  }
768  }
769  }
770  }
771  else if (geoId.subdetId() == MuonSubdetId::CSC) {
772 
773  // get station
774  CSCDetId did(geoId.rawId());
775  int station = did.station();
776  int ring = did.ring();
777 
778  //get rechits from segments by station
780  TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_begin = muonCorrelatedHitsTmp.begin();
781  TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_end = muonCorrelatedHitsTmp.end();
782 
783  for (; hits_begin!= hits_end;++hits_begin) {
784  muonCorrelatedHits.at(station-1).push_back(*hits_begin);
785  }
786 
787  if ((station == 1 && ring == 3) && dtOverlapToCheck) cscOverlapToCheck = true;
788 
789  // split 1D rechits by station
790  CSCRecHit2DCollection::range dRecHits = theCSCRecHits->get(did);
791  for (CSCRecHit2DCollection::const_iterator rechit = dRecHits.first; rechit != dRecHits.second; ++rechit) {
792 
793  if (!cscOverlapToCheck) {
794  muonRecHits.at(station-1).push_back(MuonTransientTrackingRecHit::specificBuild((&**igd),&*rechit));
795  } else {
796  tmpCSC1.push_back(MuonTransientTrackingRecHit::specificBuild((&**igd),&*rechit));
797 
798  //sort by perp, then insert to appropriate container
800  if (temp.empty()) continue;
801 
802  float center;
803  if (temp.size() > 1) {
804  center = (temp.front()->globalPosition().perp() + temp.back()->globalPosition().perp())/2.;
805  } else {
806  center = temp.front()->globalPosition().perp();
807  }
808  temp.clear();
809 
810  if (center > 550.) {
811  muonRecHits.at(2).insert(muonRecHits.at(2).end(),tmpCSC1.begin(),tmpCSC1.end());
812  } else {
813  muonRecHits.at(1).insert(muonRecHits.at(1).end(),tmpCSC1.begin(),tmpCSC1.end());
814  }
815  tmpCSC1.clear();
816  }
817  }
818  } else if (geoId.subdetId() == MuonSubdetId::RPC) {
819  LogTrace(category_) << "Wrong subdet id" << endl;
820  }
821  }//loop over GeomDets compatible with a track
822 
823  // calculate number of all and correlated hits
824  for (int stat = 0; stat < 4; stat++) {
825  theCorrelatedStationHits[stat] = muonCorrelatedHits.at(stat).size();
826  theAllStationHits[stat] = muonRecHits[stat].size();
827  }
828  LogTrace(category_) << "Hits used by the segments, by station "
829  << theCorrelatedStationHits.at(0) << " "
830  << theCorrelatedStationHits.at(1) << " "
831  << theCorrelatedStationHits.at(2) << " "
832  << theCorrelatedStationHits.at(3) << endl;
833 
834  LogTrace(category_) << "All DT 1D/CSC 2D hits, by station "
835  << theAllStationHits.at(0) << " "
836  << theAllStationHits.at(1) << " "
837  << theAllStationHits.at(2) << " "
838  << theAllStationHits.at(3) << endl;
839 
840  //station shower sizes
841  MuonTransientTrackingRecHit::MuonRecHitContainer muonRecHitsPhiTemp, muonRecHitsPhiBest;
842  TransientTrackingRecHit::ConstRecHitContainer muonRecHitsThetaTemp, muonRecHitsThetaBest;
843 
844  // send station hits to the clustering algorithm
845  for ( int stat = 0; stat != 4; stat++ ) {
846  if (!muonRecHits[stat].empty()) {
847  stable_sort(muonRecHits[stat].begin(), muonRecHits[stat].end(), LessPhi());
848 
849  float dphimax = 0;
850  for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator iseed = muonRecHits[stat].begin(); iseed != muonRecHits[stat].end(); ++iseed) {
851  if (!(*iseed)->isValid()) continue;
852  GlobalPoint refpoint = (*iseed)->globalPosition(); //starting from the one with smallest value of phi
853  muonRecHitsPhiTemp.clear();
854  muonRecHitsPhiTemp = findPhiCluster(muonRecHits[stat], refpoint); //get clustered hits for this iseed
855  if (muonRecHitsPhiTemp.size() > 1) {
856  float dphi = fabs(deltaPhi((float)muonRecHitsPhiTemp.back()->globalPosition().phi(), (float)muonRecHitsPhiTemp.front()->globalPosition().phi()));
857  if (dphi > dphimax) {
858  dphimax = dphi;
859  muonRecHitsPhiBest = muonRecHitsPhiTemp;
860  }
861  } //at least two hits
862  }//loop over seeds
863 
864  //fill showerTs
865  if (!muonRecHitsPhiBest.empty()) {
866  muonRecHits[stat] = muonRecHitsPhiBest;
867  stable_sort(muonRecHits[stat].begin(), muonRecHits[stat].end(), LessAbsMag());
868  muonRecHits[stat].front();
869  GlobalPoint refpoint = muonRecHits[stat].front()->globalPosition();
870  theStationShowerTSize.at(stat) = refpoint.mag() * dphimax;
871  }
872 
873  //for theta
874  if (!muonCorrelatedHits.at(stat).empty()) {
875 
876  float dthetamax = 0;
877  for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator iseed = muonCorrelatedHits.at(stat).begin(); iseed != muonCorrelatedHits.at(stat).end(); ++iseed) {
878  if (!(*iseed)->isValid()) continue;
879  GlobalPoint refpoint = (*iseed)->globalPosition(); //starting from the one with smallest value of phi
880  muonRecHitsThetaTemp.clear();
881  muonRecHitsThetaTemp = findThetaCluster(muonCorrelatedHits.at(stat), refpoint);
882  }//loop over seeds
883  if (muonRecHitsThetaTemp.size() > 1) {
884  float dtheta = fabs((float)muonRecHitsThetaTemp.back()->globalPosition().theta() - (float)muonRecHitsThetaTemp.front()->globalPosition().theta());
885  if (dtheta > dthetamax) {
886  dthetamax = dtheta;
887  muonRecHitsThetaBest = muonRecHitsThetaTemp;
888  }
889  } //at least two hits
890  }//not empty container2
891 
892  //fill deltaRs
893  if (muonRecHitsThetaBest.size() > 1 && muonRecHitsPhiBest.size() > 1)
894  theStationShowerDeltaR.at(stat) = sqrt(pow(muonRecHitsPhiBest.front()->globalPosition().phi()-muonRecHitsPhiBest.back()->globalPosition().phi(),2)+pow(muonRecHitsThetaBest.front()->globalPosition().theta()-muonRecHitsThetaBest.back()->globalPosition().theta(),2));
895 
896  }//not empty container
897  }//loop over station
898 
899  LogTrace(category_) << "deltaR around a track containing all the station hits, by station "
900  << theStationShowerDeltaR.at(0) << " "
901  << theStationShowerDeltaR.at(1) << " "
902  << theStationShowerDeltaR.at(2) << " "
903  << theStationShowerDeltaR.at(3) << endl;
904 
905 
906  LogTrace(category_) << "Transverse cluster size, by station "
907  << theStationShowerTSize.at(0) << " "
908  << theStationShowerTSize.at(1) << " "
909  << theStationShowerTSize.at(2) << " "
910  << theStationShowerTSize.at(3) << endl;
911 
912  return;
913 
914 }
MuonTransientTrackingRecHit::MuonRecHitContainer MuonRecHitContainer
double p() const
momentum vector magnitude
Definition: TrackBase.h:127
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
virtual void setEvent(const edm::Event &)
pass the Event to the algorithm at each event
T getParameter(std::string const &) const
int i
Definition: DBlmapReader.cc:9
std::pair< const_iterator, const_iterator > range
iterator range
Definition: RangeMap.h:50
T perp() const
Definition: PV3DBase.h:72
edm::ESHandle< MagneticField > theField
TrajectoryStateOnSurface innerStateOnSurface(const reco::Track &tk, const TrackingGeometry &geom, const MagneticField *field)
reco::MuonShower fillShowerInformation(const reco::Muon &muon, const edm::Event &, const edm::EventSetup &)
fill muon shower variables
void setServices(const edm::EventSetup &)
set the services needed
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
static const double slope[3]
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
Geom::Theta< T > theta() const
T y() const
Definition: PV3DBase.h:63
bool isGlobalMuon() const
Definition: Muon.h:218
GlobalPoint globalPosition() const
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:137
edm::EDGetTokenT< CSCSegmentCollection > theCSCSegmentsToken
edm::ESHandle< GeometricSearchTracker > theTracker
edm::ESHandle< CSCGeometry > theCSCGeometry
edm::ESHandle< DTGeometry > theDTGeometry
std::vector< int > nStationHits
number of all the muon RecHits per chamber crossed by a track (1D hits)
Definition: MuonShower.h:8
bool isStandAloneMuon() const
Definition: Muon.h:220
edm::EDGetTokenT< CSCRecHit2DCollection > theCSCRecHitToken
virtual ~MuonShowerInformationFiller()
destructor
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
static TransientTrackingRecHit::ConstRecHitContainer breakInSubRecHits(TransientTrackingRecHit::ConstRecHitPointer, int granularity)
takes a muon rechit and returns its sub-rechits given a certain granularity
std::vector< const GeomDet * > dtPositionToDets(const GlobalPoint &) const
std::vector< int > nStationCorrelatedHits
number of the muon RecHits used by segments per chamber crossed by a track
Definition: MuonShower.h:10
int iEvent
Definition: GenABIO.cc:230
static const int CSC
Definition: MuonSubdetId.h:13
T mag() const
Definition: PV3DBase.h:67
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:139
void fillHitsByStation(const reco::Muon &)
Scalar radius() const
Radius of the cylinder.
Definition: Cylinder.h:67
int n0
Definition: AMPTWrapper.h:34
static const int maxLayerId
highest layer id
Definition: DTChamberId.h:85
std::shared_ptr< TrackingRecHit const > ConstRecHitPointer
T sqrt(T t)
Definition: SSEVec.h:48
T z() const
Definition: PV3DBase.h:64
tuple result
Definition: query.py:137
virtual const BoundDisk & specificSurface() const
edm::Handle< DTRecHitCollection > theDTRecHits
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
static const int minLayerId
lowest layer id. 0 indicates a full SL
Definition: DTChamberId.h:83
DetId geographicalId() const
The label of this GeomDet.
Definition: GeomDet.h:72
GlobalPoint crossingPoint(const GlobalPoint &, const GlobalPoint &, const BarrelDetLayer *) const
#define end
Definition: vmac.h:37
static const int maxSuperLayerId
highest superlayer id
Definition: DTChamberId.h:81
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
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
double p2[4]
Definition: TauolaWrapper.h:90
#define LogTrace(id)
edm::Handle< CSCSegmentCollection > theCSCSegments
#define M_PI
virtual TrackRef outerTrack() const
reference to Track reconstructed in the muon detector only
Definition: Muon.h:51
TransientTrackingRecHit::ConstRecHitContainer findThetaCluster(TransientTrackingRecHit::ConstRecHitContainer &, const GlobalPoint &) const
edm::ESHandle< TransientTrackingRecHitBuilder > theTrackerRecHitBuilder
std::vector< ConstRecHitPointer > ConstRecHitContainer
Definition: DetId.h:18
TransientTrackingRecHit::ConstRecHitContainer hitsFromSegments(const GeomDet *, edm::Handle< DTRecSegment4DCollection >, edm::Handle< CSCSegmentCollection >) const
edm::EDGetTokenT< DTRecSegment4DCollection > theDT4DRecSegmentToken
int iseed
Definition: AMPTWrapper.h:124
std::vector< const GeomDet * > getCompatibleDets(const reco::Track &) const
edm::Handle< CSCRecHit2DCollection > theCSCRecHits
const T & get() const
Definition: EventSetup.h:55
TrajectoryStateOnSurface outerStateOnSurface(const reco::Track &tk, const TrackingGeometry &geom, const MagneticField *field)
virtual const BoundCylinder & specificSurface() const
Extension of the interface.
static const int minSuperLayerId
loweset super layer id. 0 indicates a full chamber
Definition: DTChamberId.h:79
edm::ESHandle< TransientTrackingRecHitBuilder > theMuonRecHitBuilder
std::vector< const GeomDet * > cscPositionToDets(const GlobalPoint &) const
T perp() const
Magnitude of transverse component.
#define begin
Definition: vmac.h:30
double p1[4]
Definition: TauolaWrapper.h:89
static const int RPC
Definition: MuonSubdetId.h:14
double a
Definition: hdecay.h:121
std::vector< float > stationShowerSizeT
the transverse size of the hit cluster
Definition: MuonShower.h:12
static const int DT
Definition: MuonSubdetId.h:12
edm::EDGetTokenT< DTRecHitCollection > theDTRecHitToken
edm::ESHandle< GlobalTrackingGeometry > theTrackingGeometry
Detector det() const
get the detector field from this detid
Definition: DetId.h:35
T x() const
Definition: PV3DBase.h:62
MuonRecHitContainer findPerpCluster(MuonRecHitContainer &muonRecHits) const
static MuonRecHitPointer specificBuild(const GeomDet *geom, const TrackingRecHit *rh)
std::vector< MuonRecHitPointer > MuonRecHitContainer
MuonRecHitContainer findPhiCluster(MuonRecHitContainer &, const GlobalPoint &) const
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
edm::Handle< DTRecSegment4DCollection > theDT4DRecSegments
virtual TrackRef globalTrack() const
reference to Track reconstructed in both tracked and muon detector
Definition: Muon.h:54
std::vector< float > stationShowerDeltaR
the radius of the cone containing the all the hits around the track
Definition: MuonShower.h:14
Definition: DDAxes.h:10