CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
L2MuonSeedGeneratorFromL1T.cc
Go to the documentation of this file.
1 //-------------------------------------------------
2 //
15 //
16 //--------------------------------------------------
17 
18 // Class Header
20 
21 
22 // Framework
28 
36 
39 
40 using namespace std;
41 using namespace edm;
42 using namespace l1t;
43 
44 // constructors
46  theSource(iConfig.getParameter<InputTag>("InputObjects")),
47  theL1GMTReadoutCollection(iConfig.getParameter<InputTag>("GMTReadoutCollection")), // to be removed
48  thePropagatorName(iConfig.getParameter<string>("Propagator")),
49  theL1MinPt(iConfig.getParameter<double>("L1MinPt")),
50  theL1MaxEta(iConfig.getParameter<double>("L1MaxEta")),
51  theL1MinQuality(iConfig.getParameter<unsigned int>("L1MinQuality")),
52  useOfflineSeed(iConfig.getUntrackedParameter<bool>("UseOfflineSeed", false)),
53  useUnassociatedL1(iConfig.existsAs<bool>("UseUnassociatedL1") ?
54  iConfig.getParameter<bool>("UseUnassociatedL1") : true),
55  matchingDR(iConfig.getParameter<std::vector<double> >("MatchDR")),
56  etaBins(iConfig.getParameter<std::vector<double> >("EtaMatchingBins")),
57  centralBxOnly_( iConfig.getParameter<bool>("CentralBxOnly") )
58  {
59 
60  muCollToken_ = consumes<MuonBxCollection>(theSource);
61 
62  if(useOfflineSeed) {
63  theOfflineSeedLabel = iConfig.getUntrackedParameter<InputTag>("OfflineSeedLabel");
64  offlineSeedToken_ = consumes<edm::View<TrajectorySeed> >(theOfflineSeedLabel);
65 
66  // check that number of eta bins -1 matches number of dR cones
67  if( matchingDR.size()!=etaBins.size() - 1 ) {
68  throw cms::Exception("Configuration") << "Size of MatchDR "
69  << "does not match number of eta bins." << endl;
70  }
71 
72 
73  }
74 
75  // service parameters
76  ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");
77 
78  // the services
79  theService = new MuonServiceProxy(serviceParameters);
80 
81  // the estimator
83 
84 
85  produces<L2MuonTrajectorySeedCollection>();
86 }
87 
88 // destructor
90  if (theService) delete theService;
91  if (theEstimator) delete theEstimator;
92 }
93 
94 void
97  desc.add<edm::InputTag>("GMTReadoutCollection", edm::InputTag("")); // to be removed
98  desc.add<edm::InputTag>("InputObjects",edm::InputTag("hltGmtStage2Digis"));
99  desc.add<string>("Propagator", "");
100  desc.add<double>("L1MinPt",-1.);
101  desc.add<double>("L1MaxEta",5.0);
102  desc.add<unsigned int>("L1MinQuality",0);
103  desc.addUntracked<bool>("UseOfflineSeed",false);
104  desc.add<bool>("UseUnassociatedL1", true);
105  desc.add<std::vector<double>>("MatchDR", {0.3});
106  desc.add<std::vector<double>>("EtaMatchingBins", {0., 2.5});
107  desc.add<bool>("CentralBxOnly", true);
108  desc.addUntracked<edm::InputTag>("OfflineSeedLabel", edm::InputTag(""));
109 
111  psd0.addUntracked<std::vector<std::string>>("Propagators", {
112  "SteppingHelixPropagatorAny"
113  });
114  psd0.add<bool>("RPCLayers", true);
115  psd0.addUntracked<bool>("UseMuonNavigation", true);
116  desc.add<edm::ParameterSetDescription>("ServiceParameters", psd0);
117  descriptions.add("L2MuonSeedGeneratorFromL1T",desc);
118 }
119 
121 {
122  const std::string metname = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
124 
125  auto_ptr<L2MuonTrajectorySeedCollection> output(new L2MuonTrajectorySeedCollection());
126 
128  iEvent.getByToken(muCollToken_, muColl);
129  LogTrace(metname) << "Number of muons " << muColl->size() << endl;
130 
131  edm::Handle<edm::View<TrajectorySeed> > offlineSeedHandle;
132  vector<int> offlineSeedMap;
133  if(useOfflineSeed) {
134  iEvent.getByToken(offlineSeedToken_, offlineSeedHandle);
135  LogTrace(metname) << "Number of offline seeds " << offlineSeedHandle->size() << endl;
136  offlineSeedMap = vector<int>(offlineSeedHandle->size(), 0);
137  }
138 
139  for (int ibx = muColl->getFirstBX(); ibx <= muColl->getLastBX(); ++ibx) {
140  if (centralBxOnly_ && (ibx != 0)) continue;
141  for (auto it = muColl->begin(ibx); it != muColl->end(ibx); it++){
142 
143  unsigned int quality = it->hwQual();
144  int valid_charge = it->hwChargeValid();
145 
146  float pt = it->pt();
147  float eta = it->eta();
148  float theta = 2*atan(exp(-eta));
149  float phi = it->phi();
150  int charge = it->charge();
151  // Set charge=0 for the time being if the valid charge bit is zero
152  if (!valid_charge) charge = 0;
153 
154  int link = 36 + (int)(it -> tfMuonIndex() / 3.);
155  bool barrel = true;
156  if ( (link >= 36 && link <= 41) || (link >= 66 && link <= 71)) barrel = false;
157 
158  if ( pt < theL1MinPt || fabs(eta) > theL1MaxEta ) continue;
159 
160  LogTrace(metname) << "New L2 Muon Seed" ;
161  LogTrace(metname) << "Pt = " << pt << " GeV/c";
162  LogTrace(metname) << "eta = " << eta;
163  LogTrace(metname) << "theta = " << theta << " rad";
164  LogTrace(metname) << "phi = " << phi << " rad";
165  LogTrace(metname) << "charge = " << charge;
166  LogTrace(metname) << "In Barrel? = " << barrel;
167 
168  if ( quality <= theL1MinQuality ) continue;
169  LogTrace(metname) << "quality = "<< quality;
170 
171  // Update the services
172  theService->update(iSetup);
173 
174  const DetLayer *detLayer = 0;
175  float radius = 0.;
176 
177  CLHEP::Hep3Vector vec(0.,1.,0.);
178  vec.setTheta(theta);
179  vec.setPhi(phi);
180 
181  DetId theid;
182  // Get the det layer on which the state should be put
183  if ( barrel ){
184  LogTrace(metname) << "The seed is in the barrel";
185 
186  // MB2
187  DetId id = DTChamberId(0,2,0);
188  detLayer = theService->detLayerGeometry()->idToLayer(id);
189  LogTrace(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
190 
191  const BoundSurface* sur = &(detLayer->surface());
192  const BoundCylinder* bc = dynamic_cast<const BoundCylinder*>(sur);
193 
194  radius = fabs(bc->radius()/sin(theta));
195  theid = id;
196 
197  LogTrace(metname) << "radius "<<radius;
198 
199  if ( pt < 3.5 ) pt = 3.5;
200  }
201  else {
202  LogTrace(metname) << "The seed is in the endcap";
203 
204  DetId id;
205  // ME2
206  if ( theta < Geom::pi()/2. )
207  id = CSCDetId(1,2,0,0,0);
208  else
209  id = CSCDetId(2,2,0,0,0);
210 
211  detLayer = theService->detLayerGeometry()->idToLayer(id);
212  LogTrace(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
213 
214  radius = fabs(detLayer->position().z()/cos(theta));
215  theid = id;
216 
217  if( pt < 1.0) pt = 1.0;
218  }
219 
220  vec.setMag(radius);
221 
222  GlobalPoint pos(vec.x(),vec.y(),vec.z());
223 
224  GlobalVector mom(pt*cos(phi), pt*sin(phi), pt*cos(theta)/sin(theta));
225 
226  GlobalTrajectoryParameters param(pos,mom,charge,&*theService->magneticField());
228 
229  mat[0][0] = (0.25/pt)*(0.25/pt); // sigma^2(charge/abs_momentum)
230  if ( !barrel ) mat[0][0] = (0.4/pt)*(0.4/pt);
231 
232  //Assign q/pt = 0 +- 1/pt if charge has been declared invalid
233  if (!valid_charge) mat[0][0] = (1./pt)*(1./pt);
234 
235  mat[1][1] = 0.05*0.05; // sigma^2(lambda)
236  mat[2][2] = 0.2*0.2; // sigma^2(phi)
237  mat[3][3] = 20.*20.; // sigma^2(x_transverse))
238  mat[4][4] = 20.*20.; // sigma^2(y_transverse))
239 
241 
242  const FreeTrajectoryState state(param,error);
243 
244  LogTrace(metname) << "Free trajectory State from the parameters";
245  LogTrace(metname) << debug.dumpFTS(state);
246 
247  // Propagate the state on the MB2/ME2 surface
248  TrajectoryStateOnSurface tsos = theService->propagator(thePropagatorName)->propagate(state, detLayer->surface());
249 
250  LogTrace(metname) << "State after the propagation on the layer";
251  LogTrace(metname) << debug.dumpLayer(detLayer);
252  LogTrace(metname) << debug.dumpTSOS(tsos);
253 
254  double dRcone = matchingDR[0];
255  if ( fabs(eta) < etaBins.back() ){
256  std::vector<double>::iterator lowEdge = std::upper_bound (etaBins.begin(), etaBins.end(), fabs(eta));
257  dRcone = matchingDR.at( lowEdge - etaBins.begin() - 1);
258  }
259 
260  if (tsos.isValid()) {
261 
263 
264  if(useOfflineSeed && ( !valid_charge || charge == 0) ) {
265 
266  const TrajectorySeed *assoOffseed =
267  associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, tsos, dRcone );
268 
269  if(assoOffseed!=0) {
270  PTrajectoryStateOnDet const & seedTSOS = assoOffseed->startingState();
272  tsci = assoOffseed->recHits().first,
273  tscie = assoOffseed->recHits().second;
274  for(; tsci!=tscie; ++tsci) {
275  container.push_back(*tsci);
276  }
277  output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum,
278  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) )));
279  }
280  else {
281  if(useUnassociatedL1) {
282  // convert the TSOS into a PTSOD
283  PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( tsos, theid.rawId());
284  output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum,
285  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) )));
286  }
287  }
288  }
289  else if (useOfflineSeed && valid_charge){
290  // Get the compatible dets on the layer
291  std::vector< pair<const GeomDet*,TrajectoryStateOnSurface> >
292  detsWithStates = detLayer->compatibleDets(tsos,
293  *theService->propagator(thePropagatorName),
294  *theEstimator);
295  if (detsWithStates.size()){
296 
297  TrajectoryStateOnSurface newTSOS = detsWithStates.front().second;
298  const GeomDet *newTSOSDet = detsWithStates.front().first;
299 
300  LogTrace(metname) << "Most compatible det";
301  LogTrace(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
302 
303  LogDebug(metname) << "L1 info: Det and State:";
304  LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
305 
306  if (newTSOS.isValid()){
307 
308  //LogDebug(metname) << "(x, y, z) = (" << newTSOS.globalPosition().x() << ", "
309  // << newTSOS.globalPosition().y() << ", " << newTSOS.globalPosition().z() << ")";
310  LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag() << ", phi="
311  << newTSOS.globalPosition().phi() << ", eta=" << newTSOS.globalPosition().eta() << ")";
312  LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge()*newTSOS.globalMomentum().perp() << ", phi="
313  << newTSOS.globalMomentum().phi() << ", eta=" << newTSOS.globalMomentum().eta() << ")";
314 
315  //LogDebug(metname) << "State on it";
316  //LogDebug(metname) << debug.dumpTSOS(newTSOS);
317 
318  const TrajectorySeed *assoOffseed =
319  associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, newTSOS, dRcone);
320 
321  if(assoOffseed!=0) {
322  PTrajectoryStateOnDet const & seedTSOS = assoOffseed->startingState();
324  tsci = assoOffseed->recHits().first,
325  tscie = assoOffseed->recHits().second;
326  for(; tsci!=tscie; ++tsci) {
327  container.push_back(*tsci);
328  }
329  output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum,
330  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) )));
331  }
332  else {
333  if(useUnassociatedL1) {
334  // convert the TSOS into a PTSOD
335  PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( newTSOS,newTSOSDet->geographicalId().rawId());
336  output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum,
337  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) )));
338  }
339  }
340  }
341  }
342  }
343  else {
344  // convert the TSOS into a PTSOD
345  PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( tsos, theid.rawId());
346  output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum,
347  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) )));
348  }
349  }
350  }
351  }
352 
353  iEvent.put(output);
354 }
355 
356 
357 // FIXME: does not resolve ambiguities yet!
359  std::vector<int> & offseedMap,
360  TrajectoryStateOnSurface & newTsos,
361  double dRcone ) {
362 
363  const std::string metlabel = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
364  MuonPatternRecoDumper debugtmp;
365 
366  edm::View<TrajectorySeed>::const_iterator offseed, endOffseed = offseeds->end();
367  const TrajectorySeed *selOffseed = 0;
368  double bestDr = 99999.;
369  unsigned int nOffseed(0);
370  int lastOffseed(-1);
371 
372  for(offseed=offseeds->begin(); offseed!=endOffseed; ++offseed, ++nOffseed) {
373  if(offseedMap[nOffseed]!=0) continue;
374  GlobalPoint glbPos = theService->trackingGeometry()->idToDet(offseed->startingState().detId())->surface().toGlobal(offseed->startingState().parameters().position());
375  GlobalVector glbMom = theService->trackingGeometry()->idToDet(offseed->startingState().detId())->surface().toGlobal(offseed->startingState().parameters().momentum());
376 
377  // Preliminary check
378  double preDr = deltaR( newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), glbPos.eta(), glbPos.phi() );
379  if(preDr > 1.0) continue;
380 
381  const FreeTrajectoryState offseedFTS(glbPos, glbMom, offseed->startingState().parameters().charge(), &*theService->magneticField());
382  TrajectoryStateOnSurface offseedTsos = theService->propagator(thePropagatorName)->propagate(offseedFTS, newTsos.surface());
383  LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl;
384  LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl;
385  //LogDebug(metlabel) << "(x, y, z) = (" << newTSOS.globalPosition().x() << ", "
386  // << newTSOS.globalPosition().y() << ", " << newTSOS.globalPosition().z() << ")" << std::endl;
387  LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << ", phi="
388  << offseedFTS.position().phi() << ", eta=" << offseedFTS.position().eta() << ")" << std::endl;
389  LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge()*offseedFTS.momentum().perp() << ", phi="
390  << offseedFTS.momentum().phi() << ", eta=" << offseedFTS.momentum().eta() << ")" << std::endl << std::endl;
391  //LogDebug(metlabel) << debugtmp.dumpFTS(offseedFTS) << std::endl;
392 
393  if(offseedTsos.isValid()) {
394  LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl;
395  //LogDebug(metlabel) << "(x, y, z) = (" << offseedTsos.globalPosition().x() << ", "
396  // << offseedTsos.globalPosition().y() << ", " << offseedTsos.globalPosition().z() << ")" << std::endl;
397  LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag() << ", phi="
398  << offseedTsos.globalPosition().phi() << ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl;
399  LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge()*offseedTsos.globalMomentum().perp() << ", phi="
400  << offseedTsos.globalMomentum().phi() << ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl << std::endl;
401  //LogDebug(metlabel) << debugtmp.dumpTSOS(offseedTsos) << std::endl;
402  double newDr = deltaR( newTsos.globalPosition().eta(), newTsos.globalPosition().phi(),
403  offseedTsos.globalPosition().eta(), offseedTsos.globalPosition().phi() );
404  LogDebug(metlabel) << " -- DR = " << newDr << std::endl;
405  if( newDr < dRcone && newDr<bestDr ) {
406  LogDebug(metlabel) << " --> OK! " << newDr << std::endl << std::endl;
407  selOffseed = &*offseed;
408  bestDr = newDr;
409  offseedMap[nOffseed] = 1;
410  if(lastOffseed>-1) offseedMap[lastOffseed] = 0;
411  lastOffseed = nOffseed;
412  }
413  else {
414  LogDebug(metlabel) << " --> Rejected. " << newDr << std::endl << std::endl;
415  }
416  }
417  else {
418  LogDebug(metlabel) << "Invalid offline seed TSOS after propagation!" << std::endl << std::endl;
419  }
420  }
421 
422  return selOffseed;
423 }
#define LogDebug(id)
edm::EDGetTokenT< edm::View< TrajectorySeed > > offlineSeedToken_
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
virtual const BoundSurface & surface() const =0
The surface of the GeometricSearchDet.
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
T perp() const
Definition: PV3DBase.h:72
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:462
const std::string metname
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
Geom::Theta< T > theta() const
GlobalPoint globalPosition() const
L2MuonSeedGeneratorFromL1T(const edm::ParameterSet &)
Constructor.
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
virtual std::vector< DetWithState > compatibleDets(const TrajectoryStateOnSurface &startingState, const Propagator &prop, const MeasurementEstimator &est) const
edm::Ref< MuonBxCollection > MuonRef
Definition: Muon.h:12
std::string dumpMuonId(const DetId &id) const
std::string dumpFTS(const FreeTrajectoryState &fts) const
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
std::string dumpTSOS(const TrajectoryStateOnSurface &tsos) const
void push_back(D *&d)
Definition: OwnVector.h:290
int iEvent
Definition: GenABIO.cc:230
const SurfaceType & surface() const
T mag() const
Definition: PV3DBase.h:67
const_iterator begin() const
MuonServiceProxy * theService
the event setup proxy, it takes care the services update
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:121
recHitContainer::const_iterator const_iterator
std::vector< L2MuonTrajectorySeed > L2MuonTrajectorySeedCollection
T phi() const
Definition: Phi.h:41
T z() const
Definition: PV3DBase.h:64
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
DetId geographicalId() const
The label of this GeomDet.
Definition: GeomDet.h:77
ParameterDescriptionBase * add(U const &iLabel, T const &value)
#define LogTrace(id)
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
Definition: DetId.h:18
PTrajectoryStateOnDet const & startingState() const
#define debug
Definition: HDRShower.cc:19
virtual const Surface::PositionType & position() const
Returns position of the surface.
range recHits() const
edm::EDGetTokenT< l1t::MuonBxCollection > muCollToken_
void add(std::string const &label, ParameterSetDescription const &psetDescription)
T eta() const
Definition: PV3DBase.h:76
GlobalVector globalMomentum() const
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:81
double pi()
Definition: Pi.h:31
virtual void produce(edm::Event &, const edm::EventSetup &) override
volatile std::atomic< bool > shutdown_flag false