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  centralBxOnly_( iConfig.getParameter<bool>("CentralBxOnly") )
56  {
57 
58  muCollToken_ = consumes<MuonBxCollection>(theSource);
59 
60  if(useOfflineSeed) {
61  theOfflineSeedLabel = iConfig.getUntrackedParameter<InputTag>("OfflineSeedLabel");
62  offlineSeedToken_ = consumes<edm::View<TrajectorySeed> >(theOfflineSeedLabel);
63  }
64 
65  // service parameters
66  ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");
67 
68  // the services
69  theService = new MuonServiceProxy(serviceParameters);
70 
71  // the estimator
73 
74 
75  produces<L2MuonTrajectorySeedCollection>();
76 }
77 
78 // destructor
80  if (theService) delete theService;
81  if (theEstimator) delete theEstimator;
82 }
83 
84 void
87  desc.add<edm::InputTag>("GMTReadoutCollection", edm::InputTag("")); // to be removed
88  desc.add<edm::InputTag>("InputObjects",edm::InputTag("hltGmtStage2Digis"));
89  desc.add<string>("Propagator", "");
90  desc.add<double>("L1MinPt",-1.);
91  desc.add<double>("L1MaxEta",5.0);
92  desc.add<unsigned int>("L1MinQuality",0);
93  desc.addUntracked<bool>("UseOfflineSeed",false);
94  desc.add<bool>("UseUnassociatedL1", true);
95  desc.add<bool>("CentralBxOnly", true);
96  desc.addUntracked<edm::InputTag>("OfflineSeedLabel", edm::InputTag(""));
97 
99  psd0.addUntracked<std::vector<std::string>>("Propagators", {
100  "SteppingHelixPropagatorAny"
101  });
102  psd0.add<bool>("RPCLayers", true);
103  psd0.addUntracked<bool>("UseMuonNavigation", true);
104  desc.add<edm::ParameterSetDescription>("ServiceParameters", psd0);
105  descriptions.add("L2MuonSeedGeneratorFromL1T",desc);
106 }
107 
109 {
110  const std::string metname = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
112 
113  auto_ptr<L2MuonTrajectorySeedCollection> output(new L2MuonTrajectorySeedCollection());
114 
116  iEvent.getByToken(muCollToken_, muColl);
117  LogTrace(metname) << "Number of muons " << muColl->size() << endl;
118 
119  edm::Handle<edm::View<TrajectorySeed> > offlineSeedHandle;
120  vector<int> offlineSeedMap;
121  if(useOfflineSeed) {
122  iEvent.getByToken(offlineSeedToken_, offlineSeedHandle);
123  LogTrace(metname) << "Number of offline seeds " << offlineSeedHandle->size() << endl;
124  offlineSeedMap = vector<int>(offlineSeedHandle->size(), 0);
125  }
126 
127  for (int ibx = muColl->getFirstBX(); ibx <= muColl->getLastBX(); ++ibx) {
128  if (centralBxOnly_ && (ibx != 0)) continue;
129  for (auto it = muColl->begin(ibx); it != muColl->end(ibx); it++){
130 
131  unsigned int quality = it->hwQual();
132  int valid_charge = it->hwChargeValid();
133 
134  float pt = it->pt();
135  float eta = it->eta();
136  float theta = 2*atan(exp(-eta));
137  float phi = it->phi();
138  int charge = it->charge();
139  // Set charge=0 for the time being if the valid charge bit is zero
140  if (!valid_charge) charge = 0;
141 
142  bool barrel = fabs(eta) < 1.04 ? true : false; // FIXME: to be updated once we have definition from L1
143 
144  if ( pt < theL1MinPt || fabs(eta) > theL1MaxEta ) continue;
145 
146  LogTrace(metname) << "New L2 Muon Seed" ;
147  LogTrace(metname) << "Pt = " << pt << " GeV/c";
148  LogTrace(metname) << "eta = " << eta;
149  LogTrace(metname) << "theta = " << theta << " rad";
150  LogTrace(metname) << "phi = " << phi << " rad";
151  LogTrace(metname) << "charge = " << charge;
152  LogTrace(metname) << "In Barrel? = " << barrel;
153 
154  if ( quality <= theL1MinQuality ) continue;
155  LogTrace(metname) << "quality = "<< quality;
156 
157  // Update the services
158  theService->update(iSetup);
159 
160  const DetLayer *detLayer = 0;
161  float radius = 0.;
162 
163  CLHEP::Hep3Vector vec(0.,1.,0.);
164  vec.setTheta(theta);
165  vec.setPhi(phi);
166 
167  // Get the det layer on which the state should be put
168  if ( barrel ){
169  LogTrace(metname) << "The seed is in the barrel";
170 
171  // MB2
172  DetId id = DTChamberId(0,2,0);
173  detLayer = theService->detLayerGeometry()->idToLayer(id);
174  LogTrace(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
175 
176  const BoundSurface* sur = &(detLayer->surface());
177  const BoundCylinder* bc = dynamic_cast<const BoundCylinder*>(sur);
178 
179  radius = fabs(bc->radius()/sin(theta));
180 
181  LogTrace(metname) << "radius "<<radius;
182 
183  if ( pt < 3.5 ) pt = 3.5;
184  }
185  else {
186  LogTrace(metname) << "The seed is in the endcap";
187 
188  DetId id;
189  // ME2
190  if ( theta < Geom::pi()/2. )
191  id = CSCDetId(1,2,0,0,0);
192  else
193  id = CSCDetId(2,2,0,0,0);
194 
195  detLayer = theService->detLayerGeometry()->idToLayer(id);
196  LogTrace(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
197 
198  radius = fabs(detLayer->position().z()/cos(theta));
199 
200  if( pt < 1.0) pt = 1.0;
201  }
202 
203  vec.setMag(radius);
204 
205  GlobalPoint pos(vec.x(),vec.y(),vec.z());
206 
207  GlobalVector mom(pt*cos(phi), pt*sin(phi), pt*cos(theta)/sin(theta));
208 
209  GlobalTrajectoryParameters param(pos,mom,charge,&*theService->magneticField());
211 
212  mat[0][0] = (0.25/pt)*(0.25/pt); // sigma^2(charge/abs_momentum)
213  if ( !barrel ) mat[0][0] = (0.4/pt)*(0.4/pt);
214 
215  //Assign q/pt = 0 +- 1/pt if charge has been declared invalid
216  if (!valid_charge) mat[0][0] = (1./pt)*(1./pt);
217 
218  mat[1][1] = 0.05*0.05; // sigma^2(lambda)
219  mat[2][2] = 0.2*0.2; // sigma^2(phi)
220  mat[3][3] = 20.*20.; // sigma^2(x_transverse))
221  mat[4][4] = 20.*20.; // sigma^2(y_transverse))
222 
224 
225  const FreeTrajectoryState state(param,error);
226 
227  LogTrace(metname) << "Free trajectory State from the parameters";
228  LogTrace(metname) << debug.dumpFTS(state);
229 
230  // Propagate the state on the MB2/ME2 surface
231  TrajectoryStateOnSurface tsos = theService->propagator(thePropagatorName)->propagate(state, detLayer->surface());
232 
233  LogTrace(metname) << "State after the propagation on the layer";
234  LogTrace(metname) << debug.dumpLayer(detLayer);
235  LogTrace(metname) << debug.dumpFTS(state);
236 
237  if (tsos.isValid()) {
238  // Get the compatible dets on the layer
239  std::vector< pair<const GeomDet*,TrajectoryStateOnSurface> >
240  detsWithStates = detLayer->compatibleDets(tsos,
241  *theService->propagator(thePropagatorName),
242  *theEstimator);
243  if (detsWithStates.size()){
244 
245  TrajectoryStateOnSurface newTSOS = detsWithStates.front().second;
246  const GeomDet *newTSOSDet = detsWithStates.front().first;
247 
248  LogTrace(metname) << "Most compatible det";
249  LogTrace(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
250 
251  LogDebug(metname) << "L1 info: Det and State:";
252  LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
253 
254  if (newTSOS.isValid()){
255 
256  //LogDebug(metname) << "(x, y, z) = (" << newTSOS.globalPosition().x() << ", "
257  // << newTSOS.globalPosition().y() << ", " << newTSOS.globalPosition().z() << ")";
258  LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag() << ", phi="
259  << newTSOS.globalPosition().phi() << ", eta=" << newTSOS.globalPosition().eta() << ")";
260  LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge()*newTSOS.globalMomentum().perp() << ", phi="
261  << newTSOS.globalMomentum().phi() << ", eta=" << newTSOS.globalMomentum().eta() << ")";
262 
263  //LogDebug(metname) << "State on it";
264  //LogDebug(metname) << debug.dumpTSOS(newTSOS);
265 
266  //PTrajectoryStateOnDet seedTSOS;
268 
269  if(useOfflineSeed) {
270  const TrajectorySeed *assoOffseed =
271  associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, newTSOS);
272 
273  if(assoOffseed!=0) {
274  PTrajectoryStateOnDet const & seedTSOS = assoOffseed->startingState();
276  tsci = assoOffseed->recHits().first,
277  tscie = assoOffseed->recHits().second;
278  for(; tsci!=tscie; ++tsci) {
279  container.push_back(*tsci);
280  }
281  output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum,
282  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) )));
283 
284  }
285  else {
286  if(useUnassociatedL1) {
287  // convert the TSOS into a PTSOD
288  PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( newTSOS,newTSOSDet->geographicalId().rawId());
289  output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum,
290  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) )));
291  }
292  }
293  }
294  else {
295  // convert the TSOS into a PTSOD
296  PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState( newTSOS,newTSOSDet->geographicalId().rawId());
297  output->push_back(L2MuonTrajectorySeed(seedTSOS,container,alongMomentum,
298  MuonRef(muColl, distance(muColl->begin(muColl->getFirstBX()),it) )));
299  }
300 
301  }
302  }
303  }
304  }
305 
306  }
307 
308 
309 
310  iEvent.put(output);
311 }
312 
313 
314 // FIXME: does not resolve ambiguities yet!
316  std::vector<int> & offseedMap,
317  TrajectoryStateOnSurface & newTsos) {
318 
319  const std::string metlabel = "Muon|RecoMuon|L2MuonSeedGeneratorFromL1T";
320  MuonPatternRecoDumper debugtmp;
321 
322  edm::View<TrajectorySeed>::const_iterator offseed, endOffseed = offseeds->end();
323  const TrajectorySeed *selOffseed = 0;
324  double bestDr = 99999.;
325  unsigned int nOffseed(0);
326  int lastOffseed(-1);
327 
328  for(offseed=offseeds->begin(); offseed!=endOffseed; ++offseed, ++nOffseed) {
329  if(offseedMap[nOffseed]!=0) continue;
330  GlobalPoint glbPos = theService->trackingGeometry()->idToDet(offseed->startingState().detId())->surface().toGlobal(offseed->startingState().parameters().position());
331  GlobalVector glbMom = theService->trackingGeometry()->idToDet(offseed->startingState().detId())->surface().toGlobal(offseed->startingState().parameters().momentum());
332 
333  // Preliminary check
334  double preDr = deltaR( newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), glbPos.eta(), glbPos.phi() );
335  if(preDr > 1.0) continue;
336 
337  const FreeTrajectoryState offseedFTS(glbPos, glbMom, offseed->startingState().parameters().charge(), &*theService->magneticField());
338  TrajectoryStateOnSurface offseedTsos = theService->propagator(thePropagatorName)->propagate(offseedFTS, newTsos.surface());
339  LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl;
340  LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl;
341  //LogDebug(metlabel) << "(x, y, z) = (" << newTSOS.globalPosition().x() << ", "
342  // << newTSOS.globalPosition().y() << ", " << newTSOS.globalPosition().z() << ")" << std::endl;
343  LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << ", phi="
344  << offseedFTS.position().phi() << ", eta=" << offseedFTS.position().eta() << ")" << std::endl;
345  LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge()*offseedFTS.momentum().perp() << ", phi="
346  << offseedFTS.momentum().phi() << ", eta=" << offseedFTS.momentum().eta() << ")" << std::endl << std::endl;
347  //LogDebug(metlabel) << debugtmp.dumpFTS(offseedFTS) << std::endl;
348 
349  if(offseedTsos.isValid()) {
350  LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl;
351  //LogDebug(metlabel) << "(x, y, z) = (" << offseedTsos.globalPosition().x() << ", "
352  // << offseedTsos.globalPosition().y() << ", " << offseedTsos.globalPosition().z() << ")" << std::endl;
353  LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag() << ", phi="
354  << offseedTsos.globalPosition().phi() << ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl;
355  LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge()*offseedTsos.globalMomentum().perp() << ", phi="
356  << offseedTsos.globalMomentum().phi() << ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl << std::endl;
357  //LogDebug(metlabel) << debugtmp.dumpTSOS(offseedTsos) << std::endl;
358  double newDr = deltaR( newTsos.globalPosition().eta(), newTsos.globalPosition().phi(),
359  offseedTsos.globalPosition().eta(), offseedTsos.globalPosition().phi() );
360  LogDebug(metlabel) << " -- DR = " << newDr << std::endl;
361  if( newDr<0.3 && newDr<bestDr ) { // FIXME: to be updated once we have info on eta resolution from L1
362  LogDebug(metlabel) << " --> OK! " << newDr << std::endl << std::endl;
363  selOffseed = &*offseed;
364  bestDr = newDr;
365  offseedMap[nOffseed] = 1;
366  if(lastOffseed>-1) offseedMap[lastOffseed] = 0;
367  lastOffseed = nOffseed;
368  }
369  else {
370  LogDebug(metlabel) << " --> Rejected. " << newDr << std::endl << std::endl;
371  }
372  }
373  else {
374  LogDebug(metlabel) << "Invalid offline seed TSOS after propagation!" << std::endl << std::endl;
375  }
376  }
377 
378  return selOffseed;
379 }
#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
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
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
const TrajectorySeed * associateOfflineSeedToL1(edm::Handle< edm::View< TrajectorySeed > > &, std::vector< int > &, TrajectoryStateOnSurface &)
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