CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
L2MuonSeedGeneratorFromL1T.cc
Go to the documentation of this file.
1 //-------------------------------------------------
2 //
17 //
18 //--------------------------------------------------
19 
20 // Class Header
22 
23 // Framework
30 
38 
41 
42 using namespace std;
43 using namespace edm;
44 using namespace l1t;
45 
46 //--- Functions used in output sorting
48  l1t::MuonRef Ref_L1A = A.l1tParticle();
49  l1t::MuonRef Ref_L1B = B.l1tParticle();
50  return (Ref_L1A->pt() > Ref_L1B->pt());
51 };
52 
54  l1t::MuonRef Ref_L1A = A.l1tParticle();
55  l1t::MuonRef Ref_L1B = B.l1tParticle();
56 
57  // Compare quality first
58  if (Ref_L1A->hwQual() > Ref_L1B->hwQual())
59  return true;
60  if (Ref_L1A->hwQual() < Ref_L1B->hwQual())
61  return false;
62 
63  // For same quality L1s compare pT
64  return (Ref_L1A->pt() > Ref_L1B->pt());
65 };
66 
67 // constructors
69  : theSource(iConfig.getParameter<InputTag>("InputObjects")),
70  theL1GMTReadoutCollection(iConfig.getParameter<InputTag>("GMTReadoutCollection")), // to be removed
71  thePropagatorName(iConfig.getParameter<string>("Propagator")),
72  theL1MinPt(iConfig.getParameter<double>("L1MinPt")),
73  theL1MaxEta(iConfig.getParameter<double>("L1MaxEta")),
74  theL1MinQuality(iConfig.getParameter<unsigned int>("L1MinQuality")),
75  theMinPtBarrel(iConfig.getParameter<double>("SetMinPtBarrelTo")),
76  theMinPtEndcap(iConfig.getParameter<double>("SetMinPtEndcapTo")),
77  useOfflineSeed(iConfig.getUntrackedParameter<bool>("UseOfflineSeed", false)),
78  useUnassociatedL1(iConfig.getParameter<bool>("UseUnassociatedL1")),
79  matchingDR(iConfig.getParameter<std::vector<double>>("MatchDR")),
80  etaBins(iConfig.getParameter<std::vector<double>>("EtaMatchingBins")),
81  centralBxOnly_(iConfig.getParameter<bool>("CentralBxOnly")),
82  matchType(iConfig.getParameter<unsigned int>("MatchType")),
83  sortType(iConfig.getParameter<unsigned int>("SortType")) {
84  muCollToken_ = consumes<MuonBxCollection>(theSource);
85 
86  if (useOfflineSeed) {
87  theOfflineSeedLabel = iConfig.getUntrackedParameter<InputTag>("OfflineSeedLabel");
88  offlineSeedToken_ = consumes<edm::View<TrajectorySeed>>(theOfflineSeedLabel);
89 
90  // check that number of eta bins -1 matches number of dR cones
91  if (matchingDR.size() != etaBins.size() - 1) {
92  throw cms::Exception("Configuration") << "Size of MatchDR "
93  << "does not match number of eta bins." << endl;
94  }
95  }
96 
97  if (matchType > 4 || sortType > 3) {
98  throw cms::Exception("Configuration") << "Wrong MatchType or SortType" << endl;
99  }
100 
101  // service parameters
102  ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");
103 
104  // the services
105  theService = new MuonServiceProxy(serviceParameters, consumesCollector());
106 
107  // the estimator
109 
110  produces<L2MuonTrajectorySeedCollection>();
111 }
112 
113 // destructor
115  if (theService)
116  delete theService;
117  if (theEstimator)
118  delete theEstimator;
119 }
120 
123  desc.add<edm::InputTag>("GMTReadoutCollection", edm::InputTag("")); // to be removed
124  desc.add<edm::InputTag>("InputObjects", edm::InputTag("hltGmtStage2Digis"));
125  desc.add<string>("Propagator", "");
126  desc.add<double>("L1MinPt", -1.);
127  desc.add<double>("L1MaxEta", 5.0);
128  desc.add<unsigned int>("L1MinQuality", 0);
129  desc.add<double>("SetMinPtBarrelTo", 3.5);
130  desc.add<double>("SetMinPtEndcapTo", 1.0);
131  desc.addUntracked<bool>("UseOfflineSeed", false);
132  desc.add<bool>("UseUnassociatedL1", true);
133  desc.add<std::vector<double>>("MatchDR", {0.3});
134  desc.add<std::vector<double>>("EtaMatchingBins", {0., 2.5});
135  desc.add<bool>("CentralBxOnly", true);
136  desc.add<unsigned int>("MatchType", 0)
137  ->setComment(
138  "MatchType : 0 Old matching, 1 L1 Order(1to1), 2 Min dR(1to1), 3 Higher Q(1to1), 4 All matched L1");
139  desc.add<unsigned int>("SortType", 0)->setComment("SortType : 0 not sort, 1 Pt, 2 Q and Pt");
140  desc.addUntracked<edm::InputTag>("OfflineSeedLabel", edm::InputTag(""));
141 
143  psd0.addUntracked<std::vector<std::string>>("Propagators", {"SteppingHelixPropagatorAny"});
144  psd0.add<bool>("RPCLayers", true);
145  psd0.addUntracked<bool>("UseMuonNavigation", true);
146  desc.add<edm::ParameterSetDescription>("ServiceParameters", psd0);
147  descriptions.add("L2MuonSeedGeneratorFromL1T", desc);
148 }
149 
151  const std::string metname = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
153 
154  auto output = std::make_unique<L2MuonTrajectorySeedCollection>();
155 
157  iEvent.getByToken(muCollToken_, muColl);
158  LogDebug(metname) << "Number of muons " << muColl->size() << endl;
159 
160  //--- matchType 0 : Old logic
161  if (matchType == 0) {
162  edm::Handle<edm::View<TrajectorySeed>> offlineSeedHandle;
163  vector<int> offlineSeedMap;
164  if (useOfflineSeed) {
165  iEvent.getByToken(offlineSeedToken_, offlineSeedHandle);
166  LogDebug(metname) << "Number of offline seeds " << offlineSeedHandle->size() << endl;
167  offlineSeedMap = vector<int>(offlineSeedHandle->size(), 0);
168  }
169 
170  for (int ibx = muColl->getFirstBX(); ibx <= muColl->getLastBX(); ++ibx) {
171  if (centralBxOnly_ && (ibx != 0))
172  continue;
173  for (auto it = muColl->begin(ibx); it != muColl->end(ibx); it++) {
174  unsigned int quality = it->hwQual();
175  int valid_charge = it->hwChargeValid();
176 
177  float pt = it->pt();
178  float eta = it->eta();
179  float theta = 2 * atan(exp(-eta));
180  float phi = it->phi();
181  int charge = it->charge();
182  // Set charge=0 for the time being if the valid charge bit is zero
183  if (!valid_charge)
184  charge = 0;
185 
186  // link number indices of the optical fibres that connect the uGMT with the track finders
187  //EMTF+ : 36-41, OMTF+ : 42-47, BMTF : 48-59, OMTF- : 60-65, EMTF- : 66-71
188  int link = 36 + (int)(it->tfMuonIndex() / 3.);
189  bool barrel = true;
190  if ((link >= 36 && link <= 41) || (link >= 66 && link <= 71))
191  barrel = false;
192 
193  if (pt < theL1MinPt || fabs(eta) > theL1MaxEta)
194  continue;
195 
196  LogDebug(metname) << "New L2 Muon Seed";
197  LogDebug(metname) << "Pt = " << pt << " GeV/c";
198  LogDebug(metname) << "eta = " << eta;
199  LogDebug(metname) << "theta = " << theta << " rad";
200  LogDebug(metname) << "phi = " << phi << " rad";
201  LogDebug(metname) << "charge = " << charge;
202  LogDebug(metname) << "In Barrel? = " << barrel;
203 
204  if (quality <= theL1MinQuality)
205  continue;
206  LogDebug(metname) << "quality = " << quality;
207 
208  // Update the services
209  theService->update(iSetup);
210 
211  const DetLayer *detLayer = nullptr;
212  float radius = 0.;
213 
214  CLHEP::Hep3Vector vec(0., 1., 0.);
215  vec.setTheta(theta);
216  vec.setPhi(phi);
217 
218  DetId theid;
219  // Get the det layer on which the state should be put
220  if (barrel) {
221  LogDebug(metname) << "The seed is in the barrel";
222 
223  // MB2
224  theid = DTChamberId(0, 2, 0);
225  detLayer = theService->detLayerGeometry()->idToLayer(theid);
226  LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
227 
228  const BoundSurface *sur = &(detLayer->surface());
229  const BoundCylinder *bc = dynamic_cast<const BoundCylinder *>(sur);
230 
231  radius = fabs(bc->radius() / sin(theta));
232 
233  LogDebug(metname) << "radius " << radius;
234 
235  if (pt < theMinPtBarrel)
236  pt = theMinPtBarrel;
237  } else {
238  LogDebug(metname) << "The seed is in the endcap";
239 
240  // ME2
241  theid = theta < Geom::pi() / 2. ? CSCDetId(1, 2, 0, 0, 0) : CSCDetId(2, 2, 0, 0, 0);
242 
243  detLayer = theService->detLayerGeometry()->idToLayer(theid);
244  LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
245 
246  radius = fabs(detLayer->position().z() / cos(theta));
247 
248  if (pt < theMinPtEndcap)
249  pt = theMinPtEndcap;
250  }
251 
252  vec.setMag(radius);
253 
254  GlobalPoint pos(vec.x(), vec.y(), vec.z());
255 
256  GlobalVector mom(pt * cos(phi), pt * sin(phi), pt * cos(theta) / sin(theta));
257 
258  GlobalTrajectoryParameters param(pos, mom, charge, &*theService->magneticField());
260 
261  mat[0][0] = (0.25 / pt) * (0.25 / pt); // sigma^2(charge/abs_momentum)
262  if (!barrel)
263  mat[0][0] = (0.4 / pt) * (0.4 / pt);
264 
265  //Assign q/pt = 0 +- 1/pt if charge has been declared invalid
266  if (!valid_charge)
267  mat[0][0] = (1. / pt) * (1. / pt);
268 
269  mat[1][1] = 0.05 * 0.05; // sigma^2(lambda)
270  mat[2][2] = 0.2 * 0.2; // sigma^2(phi)
271  mat[3][3] = 20. * 20.; // sigma^2(x_transverse))
272  mat[4][4] = 20. * 20.; // sigma^2(y_transverse))
273 
275 
276  const FreeTrajectoryState state(param, error);
277 
278  LogDebug(metname) << "Free trajectory State from the parameters";
279  LogDebug(metname) << debug.dumpFTS(state);
280 
281  // Propagate the state on the MB2/ME2 surface
283  theService->propagator(thePropagatorName)->propagate(state, detLayer->surface());
284 
285  LogDebug(metname) << "State after the propagation on the layer";
286  LogDebug(metname) << debug.dumpLayer(detLayer);
287  LogDebug(metname) << debug.dumpTSOS(tsos);
288 
289  double dRcone = matchingDR[0];
290  if (fabs(eta) < etaBins.back()) {
291  std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(eta));
292  dRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
293  }
294 
295  if (tsos.isValid()) {
297 
298  if (useOfflineSeed && (!valid_charge || charge == 0)) {
299  const TrajectorySeed *assoOffseed =
300  associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, tsos, dRcone);
301 
302  if (assoOffseed != nullptr) {
303  PTrajectoryStateOnDet const &seedTSOS = assoOffseed->startingState();
304  for (auto const &recHit : assoOffseed->recHits()) {
305  container.push_back(recHit);
306  }
307  output->push_back(
308  L2MuonTrajectorySeed(seedTSOS,
309  container,
311  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
312  } else {
313  if (useUnassociatedL1) {
314  // convert the TSOS into a PTSOD
316  output->push_back(
317  L2MuonTrajectorySeed(seedTSOS,
318  container,
320  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
321  }
322  }
323  } else if (useOfflineSeed && valid_charge) {
324  // Get the compatible dets on the layer
325  std::vector<pair<const GeomDet *, TrajectoryStateOnSurface>> detsWithStates =
326  detLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
327 
328  if (detsWithStates.empty() && barrel) {
329  // Fallback solution using ME2, try again to propagate but using ME2 as reference
330  DetId fallback_id;
331  theta < Geom::pi() / 2. ? fallback_id = CSCDetId(1, 2, 0, 0, 0) : fallback_id = CSCDetId(2, 2, 0, 0, 0);
332  const DetLayer *ME2DetLayer = theService->detLayerGeometry()->idToLayer(fallback_id);
333 
334  tsos = theService->propagator(thePropagatorName)->propagate(state, ME2DetLayer->surface());
335  detsWithStates =
336  ME2DetLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
337  }
338 
339  if (!detsWithStates.empty()) {
340  TrajectoryStateOnSurface newTSOS = detsWithStates.front().second;
341  const GeomDet *newTSOSDet = detsWithStates.front().first;
342 
343  LogDebug(metname) << "Most compatible det";
344  LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
345 
346  if (newTSOS.isValid()) {
347  LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag()
348  << ", phi=" << newTSOS.globalPosition().phi()
349  << ", eta=" << newTSOS.globalPosition().eta() << ")";
350  LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge() * newTSOS.globalMomentum().perp()
351  << ", phi=" << newTSOS.globalMomentum().phi()
352  << ", eta=" << newTSOS.globalMomentum().eta() << ")";
353 
354  const TrajectorySeed *assoOffseed =
355  associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, newTSOS, dRcone);
356 
357  if (assoOffseed != nullptr) {
358  PTrajectoryStateOnDet const &seedTSOS = assoOffseed->startingState();
359  for (auto const &recHit : assoOffseed->recHits()) {
360  container.push_back(recHit);
361  }
362  output->push_back(
363  L2MuonTrajectorySeed(seedTSOS,
364  container,
366  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
367  } else {
368  if (useUnassociatedL1) {
369  // convert the TSOS into a PTSOD
370  PTrajectoryStateOnDet const &seedTSOS =
372  output->push_back(
373  L2MuonTrajectorySeed(seedTSOS,
374  container,
376  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
377  }
378  }
379  }
380  }
381  } else {
382  // convert the TSOS into a PTSOD
384  output->push_back(L2MuonTrajectorySeed(
385  seedTSOS, container, alongMomentum, MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
386  }
387  }
388  }
389  }
390 
391  } // End of matchType 0
392 
393  //--- matchType > 0 : Updated logics
394  else if (matchType > 0) {
395  unsigned int nMuColl = muColl->size();
396 
397  vector<vector<double>> dRmtx;
398  vector<vector<const TrajectorySeed *>> selOffseeds;
399 
400  edm::Handle<edm::View<TrajectorySeed>> offlineSeedHandle;
401  if (useOfflineSeed) {
402  iEvent.getByToken(offlineSeedToken_, offlineSeedHandle);
403  unsigned int nOfflineSeed = offlineSeedHandle->size();
404  LogDebug(metname) << "Number of offline seeds " << nOfflineSeed << endl;
405 
406  // Initialize dRmtx and selOffseeds
407  dRmtx = vector<vector<double>>(nMuColl, vector<double>(nOfflineSeed, 999.0));
408  selOffseeds =
409  vector<vector<const TrajectorySeed *>>(nMuColl, vector<const TrajectorySeed *>(nOfflineSeed, nullptr));
410  }
411 
412  // At least one L1 muons are associated with offSeed
413  bool isAsso = false;
414 
415  //--- Fill dR matrix
416  for (int ibx = muColl->getFirstBX(); ibx <= muColl->getLastBX(); ++ibx) {
417  if (centralBxOnly_ && (ibx != 0))
418  continue;
419  for (auto it = muColl->begin(ibx); it != muColl->end(ibx); it++) {
420  //Position of given L1mu
421  unsigned int imu = distance(muColl->begin(muColl->getFirstBX()), it);
422 
423  unsigned int quality = it->hwQual();
424  int valid_charge = it->hwChargeValid();
425 
426  float pt = it->pt();
427  float eta = it->eta();
428  float theta = 2 * atan(exp(-eta));
429  float phi = it->phi();
430  int charge = it->charge();
431  // Set charge=0 for the time being if the valid charge bit is zero
432  if (!valid_charge)
433  charge = 0;
434 
435  // link number indices of the optical fibres that connect the uGMT with the track finders
436  //EMTF+ : 36-41, OMTF+ : 42-47, BMTF : 48-59, OMTF- : 60-65, EMTF- : 66-71
437  int link = 36 + (int)(it->tfMuonIndex() / 3.);
438  bool barrel = true;
439  if ((link >= 36 && link <= 41) || (link >= 66 && link <= 71))
440  barrel = false;
441 
442  if (pt < theL1MinPt || fabs(eta) > theL1MaxEta)
443  continue;
444 
445  LogDebug(metname) << "New L2 Muon Seed";
446  LogDebug(metname) << "Pt = " << pt << " GeV/c";
447  LogDebug(metname) << "eta = " << eta;
448  LogDebug(metname) << "theta = " << theta << " rad";
449  LogDebug(metname) << "phi = " << phi << " rad";
450  LogDebug(metname) << "charge = " << charge;
451  LogDebug(metname) << "In Barrel? = " << barrel;
452 
453  if (quality <= theL1MinQuality)
454  continue;
455  LogDebug(metname) << "quality = " << quality;
456 
457  // Update the services
458  theService->update(iSetup);
459 
460  const DetLayer *detLayer = nullptr;
461  float radius = 0.;
462 
463  CLHEP::Hep3Vector vec(0., 1., 0.);
464  vec.setTheta(theta);
465  vec.setPhi(phi);
466 
467  DetId theid;
468  // Get the det layer on which the state should be put
469  if (barrel) {
470  LogDebug(metname) << "The seed is in the barrel";
471 
472  // MB2
473  theid = DTChamberId(0, 2, 0);
474  detLayer = theService->detLayerGeometry()->idToLayer(theid);
475  LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
476 
477  const BoundSurface *sur = &(detLayer->surface());
478  const BoundCylinder *bc = dynamic_cast<const BoundCylinder *>(sur);
479 
480  radius = fabs(bc->radius() / sin(theta));
481 
482  LogDebug(metname) << "radius " << radius;
483 
484  if (pt < theMinPtBarrel)
485  pt = theMinPtBarrel;
486  } else {
487  LogDebug(metname) << "The seed is in the endcap";
488 
489  // ME2
490  theid = theta < Geom::pi() / 2. ? CSCDetId(1, 2, 0, 0, 0) : CSCDetId(2, 2, 0, 0, 0);
491 
492  detLayer = theService->detLayerGeometry()->idToLayer(theid);
493  LogDebug(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
494 
495  radius = fabs(detLayer->position().z() / cos(theta));
496 
497  if (pt < theMinPtEndcap)
498  pt = theMinPtEndcap;
499  }
500 
501  vec.setMag(radius);
502 
503  GlobalPoint pos(vec.x(), vec.y(), vec.z());
504 
505  GlobalVector mom(pt * cos(phi), pt * sin(phi), pt * cos(theta) / sin(theta));
506 
507  GlobalTrajectoryParameters param(pos, mom, charge, &*theService->magneticField());
509 
510  mat[0][0] = (0.25 / pt) * (0.25 / pt); // sigma^2(charge/abs_momentum)
511  if (!barrel)
512  mat[0][0] = (0.4 / pt) * (0.4 / pt);
513 
514  //Assign q/pt = 0 +- 1/pt if charge has been declared invalid
515  if (!valid_charge)
516  mat[0][0] = (1. / pt) * (1. / pt);
517 
518  mat[1][1] = 0.05 * 0.05; // sigma^2(lambda)
519  mat[2][2] = 0.2 * 0.2; // sigma^2(phi)
520  mat[3][3] = 20. * 20.; // sigma^2(x_transverse))
521  mat[4][4] = 20. * 20.; // sigma^2(y_transverse))
522 
524 
525  const FreeTrajectoryState state(param, error);
526 
527  LogDebug(metname) << "Free trajectory State from the parameters";
528  LogDebug(metname) << debug.dumpFTS(state);
529 
530  // Propagate the state on the MB2/ME2 surface
532  theService->propagator(thePropagatorName)->propagate(state, detLayer->surface());
533 
534  LogDebug(metname) << "State after the propagation on the layer";
535  LogDebug(metname) << debug.dumpLayer(detLayer);
536  LogDebug(metname) << debug.dumpTSOS(tsos);
537 
538  double dRcone = matchingDR[0];
539  if (fabs(eta) < etaBins.back()) {
540  std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(eta));
541  dRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
542  }
543 
544  if (tsos.isValid()) {
546 
547  if (useOfflineSeed) {
548  if ((!valid_charge || charge == 0)) {
549  bool isAssoOffseed = isAssociateOfflineSeedToL1(offlineSeedHandle, dRmtx, tsos, imu, selOffseeds, dRcone);
550 
551  if (isAssoOffseed) {
552  isAsso = true;
553  }
554 
555  // Using old way
556  else {
557  if (useUnassociatedL1) {
558  // convert the TSOS into a PTSOD
559  PTrajectoryStateOnDet const &seedTSOS =
561  output->push_back(
562  L2MuonTrajectorySeed(seedTSOS,
563  container,
565  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
566  }
567  }
568 
569  }
570 
571  else if (valid_charge) {
572  // Get the compatible dets on the layer
573  std::vector<pair<const GeomDet *, TrajectoryStateOnSurface>> detsWithStates =
574  detLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
575 
576  if (detsWithStates.empty() && barrel) {
577  // Fallback solution using ME2, try again to propagate but using ME2 as reference
578  DetId fallback_id;
579  theta < Geom::pi() / 2. ? fallback_id = CSCDetId(1, 2, 0, 0, 0) : fallback_id = CSCDetId(2, 2, 0, 0, 0);
580  const DetLayer *ME2DetLayer = theService->detLayerGeometry()->idToLayer(fallback_id);
581 
582  tsos = theService->propagator(thePropagatorName)->propagate(state, ME2DetLayer->surface());
583  detsWithStates =
584  ME2DetLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
585  }
586 
587  if (!detsWithStates.empty()) {
588  TrajectoryStateOnSurface newTSOS = detsWithStates.front().second;
589  const GeomDet *newTSOSDet = detsWithStates.front().first;
590 
591  LogDebug(metname) << "Most compatible det";
592  LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
593 
594  if (newTSOS.isValid()) {
595  LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag()
596  << ", phi=" << newTSOS.globalPosition().phi()
597  << ", eta=" << newTSOS.globalPosition().eta() << ")";
598  LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge() * newTSOS.globalMomentum().perp()
599  << ", phi=" << newTSOS.globalMomentum().phi()
600  << ", eta=" << newTSOS.globalMomentum().eta() << ")";
601 
602  bool isAssoOffseed =
603  isAssociateOfflineSeedToL1(offlineSeedHandle, dRmtx, newTSOS, imu, selOffseeds, dRcone);
604 
605  if (isAssoOffseed) {
606  isAsso = true;
607  }
608 
609  // Using old way
610  else {
611  if (useUnassociatedL1) {
612  // convert the TSOS into a PTSOD
613  PTrajectoryStateOnDet const &seedTSOS =
615  output->push_back(
616  L2MuonTrajectorySeed(seedTSOS,
617  container,
619  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
620  }
621  }
622  }
623  }
624  }
625  } // useOfflineSeed
626 
627  else {
628  // convert the TSOS into a PTSOD
630  output->push_back(L2MuonTrajectorySeed(
631  seedTSOS, container, alongMomentum, MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
632  }
633  }
634  }
635  }
636 
637  // MatchType : 0 Old matching, 1 L1 Order(1to1), 2 Min dR(1to1), 3 Higher Q(1to1), 4 All matched L1
638  if (useOfflineSeed && isAsso) {
639  unsigned int nOfflineSeed1 = offlineSeedHandle->size();
640  unsigned int nL1;
641  unsigned int i, j; // for the matrix element
642 
643  if (matchType == 1) {
644  vector<bool> removed_col = vector<bool>(nOfflineSeed1, false);
645 
646  for (nL1 = 0; nL1 < nMuColl; ++nL1) {
647  double tempDR = 99.;
648  unsigned int theOffs = 0;
649 
650  for (j = 0; j < nOfflineSeed1; ++j) {
651  if (removed_col[j])
652  continue;
653  if (tempDR > dRmtx[nL1][j]) {
654  tempDR = dRmtx[nL1][j];
655  theOffs = j;
656  }
657  }
658 
659  auto it = muColl->begin(muColl->getFirstBX()) + nL1;
660 
661  double newDRcone = matchingDR[0];
662  if (fabs(it->eta()) < etaBins.back()) {
663  std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
664  newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
665  }
666 
667  if (!(tempDR < newDRcone))
668  continue;
669 
670  // Remove column for given offSeed
671  removed_col[theOffs] = true;
672 
673  if (selOffseeds[nL1][theOffs] != nullptr) {
674  //put given L1mu and offSeed to output
675  edm::OwnVector<TrackingRecHit> newContainer;
676 
677  PTrajectoryStateOnDet const &seedTSOS = selOffseeds[nL1][theOffs]->startingState();
678  for (auto const &recHit : selOffseeds[nL1][theOffs]->recHits()) {
679  newContainer.push_back(recHit);
680  }
681  output->push_back(L2MuonTrajectorySeed(seedTSOS,
682  newContainer,
684  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
685  }
686  }
687  }
688 
689  else if (matchType == 2) {
690  vector<bool> removed_row = vector<bool>(nMuColl, false);
691  vector<bool> removed_col = vector<bool>(nOfflineSeed1, false);
692 
693  for (nL1 = 0; nL1 < nMuColl; ++nL1) {
694  double tempDR = 99.;
695  unsigned int theL1 = 0;
696  unsigned int theOffs = 0;
697 
698  for (i = 0; i < nMuColl; ++i) {
699  if (removed_row[i])
700  continue;
701  for (j = 0; j < nOfflineSeed1; ++j) {
702  if (removed_col[j])
703  continue;
704  if (tempDR > dRmtx[i][j]) {
705  tempDR = dRmtx[i][j];
706  theL1 = i;
707  theOffs = j;
708  }
709  }
710  }
711 
712  auto it = muColl->begin(muColl->getFirstBX()) + theL1;
713 
714  double newDRcone = matchingDR[0];
715  if (fabs(it->eta()) < etaBins.back()) {
716  std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
717  newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
718  }
719 
720  if (!(tempDR < newDRcone))
721  continue;
722 
723  // Remove row and column for given (L1mu, offSeed)
724  removed_col[theOffs] = true;
725  removed_row[theL1] = true;
726 
727  if (selOffseeds[theL1][theOffs] != nullptr) {
728  //put given L1mu and offSeed to output
729  edm::OwnVector<TrackingRecHit> newContainer;
730 
731  PTrajectoryStateOnDet const &seedTSOS = selOffseeds[theL1][theOffs]->startingState();
732  for (auto const &recHit : selOffseeds[theL1][theOffs]->recHits()) {
733  newContainer.push_back(recHit);
734  }
735  output->push_back(L2MuonTrajectorySeed(seedTSOS,
736  newContainer,
738  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
739  }
740  }
741  }
742 
743  else if (matchType == 3) {
744  vector<bool> removed_row = vector<bool>(nMuColl, false);
745  vector<bool> removed_col = vector<bool>(nOfflineSeed1, false);
746 
747  for (nL1 = 0; nL1 < nMuColl; ++nL1) {
748  double tempDR = 99.;
749  unsigned int theL1 = 0;
750  unsigned int theOffs = 0;
751  auto theit = muColl->begin(muColl->getFirstBX());
752 
753  // L1Q > 10 (L1Q = 12)
754  for (i = 0; i < nMuColl; ++i) {
755  if (removed_row[i])
756  continue;
757  theit = muColl->begin(muColl->getFirstBX()) + i;
758  if (theit->hwQual() > 10) {
759  for (j = 0; j < nOfflineSeed1; ++j) {
760  if (removed_col[j])
761  continue;
762  if (tempDR > dRmtx[i][j]) {
763  tempDR = dRmtx[i][j];
764  theL1 = i;
765  theOffs = j;
766  }
767  }
768  }
769  }
770 
771  // 6 < L1Q <= 10 (L1Q = 8)
772  if (tempDR == 99.) {
773  for (i = 0; i < nMuColl; ++i) {
774  if (removed_row[i])
775  continue;
776  theit = muColl->begin(muColl->getFirstBX()) + i;
777  if ((theit->hwQual() <= 10) && (theit->hwQual() > 6)) {
778  for (j = 0; j < nOfflineSeed1; ++j) {
779  if (removed_col[j])
780  continue;
781  if (tempDR > dRmtx[i][j]) {
782  tempDR = dRmtx[i][j];
783  theL1 = i;
784  theOffs = j;
785  }
786  }
787  }
788  }
789  }
790 
791  // L1Q <= 6 (L1Q = 4)
792  if (tempDR == 99.) {
793  for (i = 0; i < nMuColl; ++i) {
794  if (removed_row[i])
795  continue;
796  theit = muColl->begin(muColl->getFirstBX()) + i;
797  if (theit->hwQual() <= 6) {
798  for (j = 0; j < nOfflineSeed1; ++j) {
799  if (removed_col[j])
800  continue;
801  if (tempDR > dRmtx[i][j]) {
802  tempDR = dRmtx[i][j];
803  theL1 = i;
804  theOffs = j;
805  }
806  }
807  }
808  }
809  }
810 
811  auto it = muColl->begin(muColl->getFirstBX()) + theL1;
812 
813  double newDRcone = matchingDR[0];
814  if (fabs(it->eta()) < etaBins.back()) {
815  std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
816  newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
817  }
818 
819  if (!(tempDR < newDRcone))
820  continue;
821 
822  // Remove row and column for given (L1mu, offSeed)
823  removed_col[theOffs] = true;
824  removed_row[theL1] = true;
825 
826  if (selOffseeds[theL1][theOffs] != nullptr) {
827  //put given L1mu and offSeed to output
828  edm::OwnVector<TrackingRecHit> newContainer;
829 
830  PTrajectoryStateOnDet const &seedTSOS = selOffseeds[theL1][theOffs]->startingState();
831  for (auto const &recHit : selOffseeds[theL1][theOffs]->recHits()) {
832  newContainer.push_back(recHit);
833  }
834  output->push_back(L2MuonTrajectorySeed(seedTSOS,
835  newContainer,
837  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
838  }
839  }
840  }
841 
842  else if (matchType == 4) {
843  for (i = 0; i < nMuColl; ++i) {
844  auto it = muColl->begin(muColl->getFirstBX()) + i;
845 
846  double tempDR = 99.;
847  unsigned int theOffs = 0;
848 
849  double newDRcone = matchingDR[0];
850  if (fabs(it->eta()) < etaBins.back()) {
851  std::vector<double>::iterator lowEdge = std::upper_bound(etaBins.begin(), etaBins.end(), fabs(it->eta()));
852  newDRcone = matchingDR.at(lowEdge - etaBins.begin() - 1);
853  }
854 
855  for (j = 0; j < nOfflineSeed1; ++j) {
856  if (tempDR > dRmtx[i][j]) {
857  tempDR = dRmtx[i][j];
858  theOffs = j;
859  }
860  }
861 
862  if (!(tempDR < newDRcone))
863  continue;
864 
865  if (selOffseeds[i][theOffs] != nullptr) {
866  //put given L1mu and offSeed to output
867  edm::OwnVector<TrackingRecHit> newContainer;
868 
869  PTrajectoryStateOnDet const &seedTSOS = selOffseeds[i][theOffs]->startingState();
870  for (auto const &recHit : selOffseeds[i][theOffs]->recHits()) {
871  newContainer.push_back(recHit);
872  }
873  output->push_back(L2MuonTrajectorySeed(seedTSOS,
874  newContainer,
876  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()), it))));
877  }
878  }
879  }
880 
881  } // useOfflineSeed && isAsso
882 
883  } // matchType > 0
884 
885  //--- SortType 1 : by L1 pT
886  if (sortType == 1) {
887  sort(output->begin(), output->end(), sortByL1Pt);
888  }
889 
890  //--- SortType 2 : by L1 quality and pT
891  else if (sortType == 2) {
892  sort(output->begin(), output->end(), sortByL1QandPt);
893  }
894 
895  iEvent.put(std::move(output));
896 }
897 
898 // FIXME: does not resolve ambiguities yet!
901  std::vector<int> &offseedMap,
902  TrajectoryStateOnSurface &newTsos,
903  double dRcone) {
904  const std::string metlabel = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
905  MuonPatternRecoDumper debugtmp;
906 
907  edm::View<TrajectorySeed>::const_iterator offseed, endOffseed = offseeds->end();
908  const TrajectorySeed *selOffseed = nullptr;
909  double bestDr = 99999.;
910  unsigned int nOffseed(0);
911  int lastOffseed(-1);
912 
913  for (offseed = offseeds->begin(); offseed != endOffseed; ++offseed, ++nOffseed) {
914  if (offseedMap[nOffseed] != 0)
915  continue;
916  GlobalPoint glbPos = theService->trackingGeometry()
917  ->idToDet(offseed->startingState().detId())
918  ->surface()
919  .toGlobal(offseed->startingState().parameters().position());
920  GlobalVector glbMom = theService->trackingGeometry()
921  ->idToDet(offseed->startingState().detId())
922  ->surface()
923  .toGlobal(offseed->startingState().parameters().momentum());
924 
925  // Preliminary check
926  double preDr = deltaR(newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), glbPos.eta(), glbPos.phi());
927  if (preDr > 1.0)
928  continue;
929 
930  const FreeTrajectoryState offseedFTS(
931  glbPos, glbMom, offseed->startingState().parameters().charge(), &*theService->magneticField());
932  TrajectoryStateOnSurface offseedTsos =
933  theService->propagator(thePropagatorName)->propagate(offseedFTS, newTsos.surface());
934  LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl;
935  LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl;
936  LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << ", phi=" << offseedFTS.position().phi()
937  << ", eta=" << offseedFTS.position().eta() << ")" << std::endl;
938  LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge() * offseedFTS.momentum().perp()
939  << ", phi=" << offseedFTS.momentum().phi() << ", eta=" << offseedFTS.momentum().eta() << ")"
940  << std::endl
941  << std::endl;
942 
943  if (offseedTsos.isValid()) {
944  LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl;
945  LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag()
946  << ", phi=" << offseedTsos.globalPosition().phi()
947  << ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl;
948  LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge() * offseedTsos.globalMomentum().perp()
949  << ", phi=" << offseedTsos.globalMomentum().phi()
950  << ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl
951  << std::endl;
952  double newDr = deltaR(newTsos.globalPosition().eta(),
953  newTsos.globalPosition().phi(),
954  offseedTsos.globalPosition().eta(),
955  offseedTsos.globalPosition().phi());
956  LogDebug(metlabel) << " -- DR = " << newDr << std::endl;
957  if (newDr < dRcone && newDr < bestDr) {
958  LogDebug(metlabel) << " --> OK! " << newDr << std::endl << std::endl;
959  selOffseed = &*offseed;
960  bestDr = newDr;
961  offseedMap[nOffseed] = 1;
962  if (lastOffseed > -1)
963  offseedMap[lastOffseed] = 0;
964  lastOffseed = nOffseed;
965  } else {
966  LogDebug(metlabel) << " --> Rejected. " << newDr << std::endl << std::endl;
967  }
968  } else {
969  LogDebug(metlabel) << "Invalid offline seed TSOS after propagation!" << std::endl << std::endl;
970  }
971  }
972 
973  return selOffseed;
974 }
975 
978  std::vector<std::vector<double>> &dRmtx,
979  TrajectoryStateOnSurface &newTsos,
980  unsigned int imu,
981  std::vector<std::vector<const TrajectorySeed *>> &selOffseeds,
982  double dRcone) {
983  const std::string metlabel = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
984  bool isAssociated = false;
985  MuonPatternRecoDumper debugtmp;
986 
987  edm::View<TrajectorySeed>::const_iterator offseed, endOffseed = offseeds->end();
988  unsigned int nOffseed(0);
989 
990  for (offseed = offseeds->begin(); offseed != endOffseed; ++offseed, ++nOffseed) {
991  GlobalPoint glbPos = theService->trackingGeometry()
992  ->idToDet(offseed->startingState().detId())
993  ->surface()
994  .toGlobal(offseed->startingState().parameters().position());
995  GlobalVector glbMom = theService->trackingGeometry()
996  ->idToDet(offseed->startingState().detId())
997  ->surface()
998  .toGlobal(offseed->startingState().parameters().momentum());
999 
1000  // Preliminary check
1001  double preDr = deltaR(newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), glbPos.eta(), glbPos.phi());
1002  if (preDr > 1.0)
1003  continue;
1004 
1005  const FreeTrajectoryState offseedFTS(
1006  glbPos, glbMom, offseed->startingState().parameters().charge(), &*theService->magneticField());
1007  TrajectoryStateOnSurface offseedTsos =
1008  theService->propagator(thePropagatorName)->propagate(offseedFTS, newTsos.surface());
1009  LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl;
1010  LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl;
1011  LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << ", phi=" << offseedFTS.position().phi()
1012  << ", eta=" << offseedFTS.position().eta() << ")" << std::endl;
1013  LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge() * offseedFTS.momentum().perp()
1014  << ", phi=" << offseedFTS.momentum().phi() << ", eta=" << offseedFTS.momentum().eta() << ")"
1015  << std::endl
1016  << std::endl;
1017 
1018  if (offseedTsos.isValid()) {
1019  LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl;
1020  LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag()
1021  << ", phi=" << offseedTsos.globalPosition().phi()
1022  << ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl;
1023  LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge() * offseedTsos.globalMomentum().perp()
1024  << ", phi=" << offseedTsos.globalMomentum().phi()
1025  << ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl
1026  << std::endl;
1027  double newDr = deltaR(newTsos.globalPosition().eta(),
1028  newTsos.globalPosition().phi(),
1029  offseedTsos.globalPosition().eta(),
1030  offseedTsos.globalPosition().phi());
1031 
1032  LogDebug(metlabel) << " -- DR = " << newDr << std::endl;
1033  if (newDr < dRcone) {
1034  LogDebug(metlabel) << " --> OK! " << newDr << std::endl << std::endl;
1035 
1036  dRmtx[imu][nOffseed] = newDr;
1037  selOffseeds[imu][nOffseed] = &*offseed;
1038 
1039  isAssociated = true;
1040  } else {
1041  LogDebug(metlabel) << " --> Rejected. " << newDr << std::endl << std::endl;
1042  }
1043  } else {
1044  LogDebug(metlabel) << "Invalid offline seed TSOS after propagation!" << std::endl << std::endl;
1045  }
1046  }
1047 
1048  return isAssociated;
1049 }
edm::EDGetTokenT< edm::View< TrajectorySeed > > offlineSeedToken_
T getUntrackedParameter(std::string const &, T const &) const
virtual const BoundSurface & surface() const =0
The surface of the GeometricSearchDet.
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:133
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
T perp() const
Definition: PV3DBase.h:69
ParameterDescriptionBase * addUntracked(U const &iLabel, T const &value)
bool centralBxOnly_
use central bx only muons
std::string dumpLayer(const DetLayer *layer) const
const TrajectorySeed * associateOfflineSeedToL1(edm::Handle< edm::View< TrajectorySeed > > &, std::vector< int > &, TrajectoryStateOnSurface &, double)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:539
bool sortByL1Pt(L2MuonTrajectorySeed &A, L2MuonTrajectorySeed &B)
const std::string metname
uint32_t const *__restrict__ Quality * quality
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
Geom::Theta< T > theta() const
l1t::MuonRef l1tParticle() const
__host__ __device__ constexpr RandomIt upper_bound(RandomIt first, RandomIt last, const T &value, Compare comp={})
GlobalPoint globalPosition() const
L2MuonSeedGeneratorFromL1T(const edm::ParameterSet &)
Constructor.
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
T1 phi() const
Definition: Phi.h:78
virtual std::vector< DetWithState > compatibleDets(const TrajectoryStateOnSurface &startingState, const Propagator &prop, const MeasurementEstimator &est) const
std::string dumpMuonId(const DetId &id) const
std::string dumpFTS(const FreeTrajectoryState &fts) const
~L2MuonSeedGeneratorFromL1T() override
Destructor.
std::string dumpTSOS(const TrajectoryStateOnSurface &tsos) const
void push_back(D *&d)
Definition: OwnVector.h:326
int iEvent
Definition: GenABIO.cc:224
const SurfaceType & surface() const
T mag() const
Definition: PV3DBase.h:64
const_iterator begin() const
MuonServiceProxy * theService
the event setup proxy, it takes care the services update
T z() const
Definition: PV3DBase.h:61
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
bool isAssociateOfflineSeedToL1(edm::Handle< edm::View< TrajectorySeed > > &, std::vector< std::vector< double > > &, TrajectoryStateOnSurface &, unsigned int, std::vector< std::vector< const TrajectorySeed * > > &, double)
def move
Definition: eostools.py:511
DetId geographicalId() const
The label of this GeomDet.
Definition: GeomDet.h:64
bool sortByL1QandPt(L2MuonTrajectorySeed &A, L2MuonTrajectorySeed &B)
static const std::string B
ParameterDescriptionBase * add(U const &iLabel, T const &value)
RecHitRange recHits() const
Definition: DetId.h:17
PTrajectoryStateOnDet const & startingState() const
#define debug
Definition: HDRShower.cc:19
virtual const Surface::PositionType & position() const
Returns position of the surface.
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
edm::EDGetTokenT< l1t::MuonBxCollection > muCollToken_
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
thePropagatorName(iConfig.getParameter< std::string >("Propagator"))
void add(std::string const &label, ParameterSetDescription const &psetDescription)
T eta() const
Definition: PV3DBase.h:73
GlobalVector globalMomentum() const
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:86
void produce(edm::Event &, const edm::EventSetup &) override
edm::Ref< MuonBxCollection > MuonRef
Definition: Muon.h:13
constexpr double pi()
Definition: Pi.h:31
#define LogDebug(id)