CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DTHitAssociator.cc
Go to the documentation of this file.
11 #include <string>
12 #include <sstream>
13 
14 using namespace std;
15 
16 // Constructor
18  edm::ConsumesCollector && iC) :
19  DTsimhitsTag(conf.getParameter<edm::InputTag>("DTsimhitsTag")),
20  DTsimhitsXFTag(conf.getParameter<edm::InputTag>("DTsimhitsXFTag")),
21  DTdigiTag(conf.getParameter<edm::InputTag>("DTdigiTag")),
22  DTdigisimlinkTag(conf.getParameter<edm::InputTag>("DTdigisimlinkTag")),
23  DTrechitTag(conf.getParameter<edm::InputTag>("DTrechitTag")),
24 
25  // nice printout of DT hits
26  dumpDT(conf.getParameter<bool>("dumpDT")),
27  // CrossingFrame used or not ?
28  crossingframe(conf.getParameter<bool>("crossingframe")),
29  // Event contain the DTDigiSimLink collection ?
30  links_exist(conf.getParameter<bool>("links_exist")),
31  // associatorByWire links to a RecHit all the "valid" SimHits on the same DT wire
32  associatorByWire(conf.getParameter<bool>("associatorByWire")),
33 
34  printRtS(true)
35 {
36 
37  if ( crossingframe) {
39  }
40  else if (!DTsimhitsTag.label().empty()) {
41  iC.consumes<edm::PSimHitContainer>(DTsimhitsTag);
42  }
43  iC.consumes<DTDigiCollection>(DTdigiTag);
45 
46  if ( dumpDT && printRtS ) {
47  iC.consumes<DTRecHitCollection>(DTrechitTag);
48  }
49 
50 }
51 
52 
54  // input collection labels
55  DTsimhitsTag(conf.getParameter<edm::InputTag>("DTsimhitsTag")),
56  DTsimhitsXFTag(conf.getParameter<edm::InputTag>("DTsimhitsXFTag")),
57  DTdigiTag(conf.getParameter<edm::InputTag>("DTdigiTag")),
58  DTdigisimlinkTag(conf.getParameter<edm::InputTag>("DTdigisimlinkTag")),
59  DTrechitTag(conf.getParameter<edm::InputTag>("DTrechitTag")),
60 
61  // nice printout of DT hits
62  dumpDT(conf.getParameter<bool>("dumpDT")),
63  // CrossingFrame used or not ?
64  crossingframe(conf.getParameter<bool>("crossingframe")),
65  // Event contain the DTDigiSimLink collection ?
66  links_exist(conf.getParameter<bool>("links_exist")),
67  // associatorByWire links to a RecHit all the "valid" SimHits on the same DT wire
68  associatorByWire(conf.getParameter<bool>("associatorByWire")),
69 
70  printRtS(true)
71 
72 {
73  initEvent(iEvent,iSetup);
74 }
75 
77 
78  LogTrace("DTHitAssociator") <<"DTHitAssociator constructor: dumpDT = "<<dumpDT
79  <<", crossingframe = "<<crossingframe<<", links_exist = "<<links_exist
80  <<", associatorByWire = "<<associatorByWire;
81 
82  if (!links_exist && !associatorByWire) {
83  edm::LogWarning("DTHitAssociator")<<"*** WARNING in DTHitAssociator::DTHitAssociator: associatorByWire reset to TRUE !"
84  <<" \t (missing DTDigiSimLinkCollection ?)";
85  associatorByWire = true;
86  }
87 
88  // need DT Geometry to discard hits for which the drift time parametrisation is not applicable
90  iSetup.get<MuonGeometryRecord>().get(muonGeom);
91 
92  // Get the DT SimHits from the event and map PSimHit by DTWireId
93  mapOfSimHit.clear();
94  bool takeHit(true);
95 
96  if (crossingframe) {
98  LogTrace("DTHitAssociator") <<"getting CrossingFrame<PSimHit> collection - "<<DTsimhitsXFTag;
99  iEvent.getByLabel(DTsimhitsXFTag,xFrame);
100  auto_ptr<MixCollection<PSimHit> >
101  DTsimhits( new MixCollection<PSimHit>(xFrame.product()) );
102  LogTrace("DTHitAssociator") <<"... size = "<<DTsimhits->size();
104  for (isimhit=DTsimhits->begin(); isimhit!= DTsimhits->end(); isimhit++) {
105  DTWireId wireid((*isimhit).detUnitId());
106  takeHit = SimHitOK(muonGeom, *isimhit);
107  mapOfSimHit[wireid].push_back(make_pair(*isimhit,takeHit));
108  }
109  }
110  else if (!DTsimhitsTag.label().empty()) {
112  LogTrace("DTHitAssociator") <<"getting PSimHit collection - "<<DTsimhitsTag;
113  iEvent.getByLabel(DTsimhitsTag,DTsimhits);
114  LogTrace("DTHitAssociator") <<"... size = "<<DTsimhits->size();
115  edm::PSimHitContainer::const_iterator isimhit;
116  for (isimhit=DTsimhits->begin(); isimhit!= DTsimhits->end(); isimhit++) {
117  DTWireId wireid((*isimhit).detUnitId());
118  takeHit = SimHitOK(muonGeom, *isimhit);
119  mapOfSimHit[wireid].push_back(make_pair(*isimhit,takeHit));
120  }
121  }
122 
123  // Get the DT Digi collection from the event
124  mapOfDigi.clear();
126  LogTrace("DTHitAssociator") <<"getting DTDigi collection - "<<DTdigiTag;
127  iEvent.getByLabel(DTdigiTag,digis);
128 
129  if (digis.isValid()) {
130  // Map DTDigi by DTWireId
131  for (DTDigiCollection::DigiRangeIterator detUnit=digis->begin(); detUnit !=digis->end(); ++detUnit) {
132  const DTLayerId& layerid = (*detUnit).first;
133  const DTDigiCollection::Range& range = (*detUnit).second;
134 
136  for (digi = range.first; digi != range.second; ++digi){
137  DTWireId wireid(layerid,(*digi).wire());
138  mapOfDigi[wireid].push_back(*digi);
139  }
140  }
141  } else {
142  LogTrace("DTHitAssociator") <<"... NO DTDigi collection found !";
143  }
144 
145  mapOfLinks.clear();
146  if (links_exist) {
147  // Get the DT DigiSimLink collection from the event and map DTDigiSimLink by DTWireId
149  LogTrace("DTHitAssociator") <<"getting DTDigiSimLink collection - "<<DTdigisimlinkTag;
150  iEvent.getByLabel(DTdigisimlinkTag,digisimlinks);
151 
152  for (DTDigiSimLinkCollection::DigiRangeIterator detUnit=digisimlinks->begin();
153  detUnit !=digisimlinks->end();
154  ++detUnit) {
155  const DTLayerId& layerid = (*detUnit).first;
156  const DTDigiSimLinkCollection::Range& range = (*detUnit).second;
157 
159  for (link=range.first; link!=range.second; ++link){
160  DTWireId wireid(layerid,(*link).wire());
161  mapOfLinks[wireid].push_back(*link);
162  }
163  }
164  }
165 
166  if(dumpDT && printRtS) {
167 
168  // Get the DT rechits from the event
170  LogTrace("DTHitAssociator") <<"getting DTRecHit1DPair collection - "<<DTrechitTag;
171  iEvent.getByLabel(DTrechitTag,DTrechits);
172  LogTrace("DTHitAssociator") <<"... size = "<<DTrechits->size();
173 
174  // map DTRecHit1DPair by DTWireId
175  mapOfRecHit.clear();
177  for(rechit=DTrechits->begin(); rechit!=DTrechits->end(); ++rechit) {
178  DTWireId wireid = (*rechit).wireId();
179  mapOfRecHit[wireid].push_back(*rechit);
180  }
181 
182  if (mapOfSimHit.end() != mapOfSimHit.begin()) {
183  edm::LogVerbatim("DTHitAssociator")<<"\n *** Dump DT PSimHit's ***";
184 
185  int jwire = 0;
186  int ihit = 0;
187 
188  for(SimHitMap::const_iterator mapIT=mapOfSimHit.begin();
189  mapIT!=mapOfSimHit.end();
190  ++mapIT , jwire++) {
191 
192  DTWireId wireid = (*mapIT).first;
193  for (vector<PSimHit_withFlag>::const_iterator hitIT = mapOfSimHit[wireid].begin();
194  hitIT != mapOfSimHit[wireid].end();
195  hitIT++ , ihit++) {
196  PSimHit hit = hitIT->first;
197  edm::LogVerbatim("DTHitAssociator")
198  <<"PSimHit "<<ihit <<", wire "<<wireid //<<", detID = "<<hit.detUnitId()
199  <<", SimTrack Id:"<<hit.trackId()<<"/Evt:(" <<hit.eventId().event()<<","<<hit.eventId().bunchCrossing()<<") "
200  <<", pdg = "<<hit.particleType()<<", procs = "<<hit.processType();
201  }
202  }
203  } else {
204  edm::LogVerbatim("DTHitAssociator")<<"\n *** There are NO DT PSimHit's ***";
205  }
206 
207  if(mapOfDigi.end() != mapOfDigi.begin()) {
208 
209  int jwire = 0;
210  int ihit = 0;
211 
212  for(DigiMap::const_iterator mapIT=mapOfDigi.begin(); mapIT!=mapOfDigi.end(); ++mapIT , jwire++) {
213  edm::LogVerbatim("DTHitAssociator")<<"\n *** Dump DT digis ***";
214 
215  DTWireId wireid = (*mapIT).first;
216  for (vector<DTDigi>::const_iterator hitIT = mapOfDigi[wireid].begin();
217  hitIT != mapOfDigi[wireid].end();
218  hitIT++ , ihit++) {
219  edm::LogVerbatim("DTHitAssociator")
220  <<"DTDigi "<<ihit<<", wire "<<wireid<<", number = "<<hitIT->number()<<", TDC counts = "<<hitIT->countsTDC();
221  }
222  }
223  } else {
224  LogTrace("DTHitAssociator")<<"\n *** There are NO DTDigi's ***";
225  }
226 
227  if (mapOfRecHit.end() != mapOfRecHit.begin())
228  edm::LogVerbatim("DTHitAssociator")<<"\n *** Analyze DTRecHitCollection by DTWireId ***";
229 
230  int iwire = 0;
231  for(RecHitMap::const_iterator mapIT=mapOfRecHit.begin();
232  mapIT!=mapOfRecHit.end();
233  ++mapIT , iwire++) {
234 
235  DTWireId wireid = (*mapIT).first;
236  edm::LogVerbatim("DTHitAssociator")<<"\n===================================================================";
237  edm::LogVerbatim("DTHitAssociator")<<"wire index = "<<iwire<<" *** DTWireId = "<<" ("<<wireid<<")";
238 
239  if(mapOfSimHit.find(wireid) != mapOfSimHit.end()) {
240  edm::LogVerbatim("DTHitAssociator")<<"\n"<<mapOfSimHit[wireid].size()<<" SimHits (PSimHit):";
241 
242  for (vector<PSimHit_withFlag>::const_iterator hitIT = mapOfSimHit[wireid].begin();
243  hitIT != mapOfSimHit[wireid].end();
244  ++hitIT) {
245  stringstream tId;
246  tId << (hitIT->first).trackId();
247  string simhitlog = "\t SimTrack Id = "+tId.str();
248  if (hitIT->second) simhitlog = simhitlog + "\t -VALID HIT-";
249  else simhitlog = simhitlog + "\t -not valid hit-";
250  edm::LogVerbatim("DTHitAssociator")<<simhitlog;
251  }
252  }
253 
254  if(mapOfLinks.find(wireid) != mapOfLinks.end()) {
255  edm::LogVerbatim("DTHitAssociator")<<"\n"<<mapOfLinks[wireid].size()<<" Links (DTDigiSimLink):";
256 
257  for (vector<DTDigiSimLink>::const_iterator hitIT = mapOfLinks[wireid].begin();
258  hitIT != mapOfLinks[wireid].end();
259  ++hitIT) {
260  edm::LogVerbatim("DTHitAssociator")
261  <<"\t digi number = "<<hitIT->number()<<", time = "<<hitIT->time()<<", SimTrackId = "<<hitIT->SimTrackId();
262  }
263  }
264 
265  if(mapOfDigi.find(wireid) != mapOfDigi.end()) {
266  edm::LogVerbatim("DTHitAssociator")<<"\n"<<mapOfDigi[wireid].size()<<" Digis (DTDigi):";
267  for (vector<DTDigi>::const_iterator hitIT = mapOfDigi[wireid].begin();
268  hitIT != mapOfDigi[wireid].end();
269  ++hitIT) {
270  edm::LogVerbatim("DTHitAssociator")<<"\t digi number = "<<hitIT->number()<<", time = "<<hitIT->time();
271  }
272  }
273 
274  edm::LogVerbatim("DTHitAssociator")<<"\n"<<(*mapIT).second.size()<<" RecHits (DTRecHit1DPair):";
275  for(vector<DTRecHit1DPair>::const_iterator vIT =(*mapIT).second.begin();
276  vIT !=(*mapIT).second.end();
277  ++vIT) {
278  edm::LogVerbatim("DTHitAssociator")<<"\t digi time = "<<vIT->digiTime();
279  }
280  }
281  }
282 }
283 // end of constructor
284 
285 std::vector<DTHitAssociator::SimHitIdpr> DTHitAssociator::associateHitId(const TrackingRecHit & hit) const {
286 
287  std::vector<SimHitIdpr> simtrackids;
288  const TrackingRecHit * hitp = &hit;
289  const DTRecHit1D * dtrechit = dynamic_cast<const DTRecHit1D *>(hitp);
290 
291  if (dtrechit) {
292  simtrackids = associateDTHitId(dtrechit);
293  } else {
294  edm::LogWarning("DTHitAssociator")<<"*** WARNING in DTHitAssociator::associateHitId, null dynamic_cast !";
295  }
296  return simtrackids;
297 }
298 
299 std::vector<DTHitAssociator::SimHitIdpr> DTHitAssociator::associateDTHitId(const DTRecHit1D * dtrechit) const {
300 
301  std::vector<SimHitIdpr> matched;
302 
303  DTWireId wireid = dtrechit->wireId();
304 
305  if(associatorByWire) {
306  // matching based on DTWireId : take only "valid" SimHits on that wire
307 
308  auto found = mapOfSimHit.find(wireid);
309  if(found != mapOfSimHit.end()) {
310  for (vector<PSimHit_withFlag>::const_iterator hitIT = found->second.begin();
311  hitIT != found->second.end();
312  ++hitIT) {
313 
314  bool valid_hit = hitIT->second;
315  PSimHit this_hit = hitIT->first;
316 
317  if (valid_hit) {
318  SimHitIdpr currentId(this_hit.trackId(), this_hit.eventId());
319  matched.push_back(currentId);
320  }
321  }
322  }
323  }
324 
325  else {
326  // matching based on DTDigiSimLink
327 
328  float theTime = dtrechit->digiTime();
329  int theNumber(-1);
330 
331  auto found = mapOfLinks.find(wireid);
332 
333  if (found != mapOfLinks.end()) {
334  // DTDigiSimLink::time() is set equal to DTDigi::time() only for the DTDigiSimLink corresponding to the digitizied PSimHit
335  // other DTDigiSimLinks associated to the same DTDigi are identified by having the same DTDigiSimLink::number()
336 
337  // first find the associated digi Number
338  for (vector<DTDigiSimLink>::const_iterator linkIT = found->second.begin();
339  linkIT != found->second.end();
340  ++linkIT ) {
341  float digitime = linkIT->time();
342  if (fabs(digitime-theTime)<0.1) {
343  theNumber = linkIT->number();
344  }
345  }
346 
347  // then get all the DTDigiSimLinks with that digi Number (corresponding to valid SimHits
348  // within a time window of the order of the time resolution, specified in the DTDigitizer)
349  for (vector<DTDigiSimLink>::const_iterator linkIT = found->second.begin();
350  linkIT != found->second.end();
351  ++linkIT ) {
352 
353  int digiNr = linkIT->number();
354  if (digiNr == theNumber) {
355  SimHitIdpr currentId(linkIT->SimTrackId(), linkIT->eventId());
356  matched.push_back(currentId);
357  }
358  }
359 
360  } else {
361  edm::LogError("DTHitAssociator")<<"*** ERROR in DTHitAssociator::associateDTHitId - DTRecHit1D: "
362  <<*dtrechit<<" has no associated DTDigiSimLink !"<<endl;
363  return matched;
364  }
365  }
366 
367  return matched;
368 }
369 
370 std::vector<PSimHit> DTHitAssociator::associateHit(const TrackingRecHit & hit) const {
371 
372  std::vector<PSimHit> simhits;
373  std::vector<SimHitIdpr> simtrackids;
374 
375  const TrackingRecHit * hitp = &hit;
376  const DTRecHit1D * dtrechit = dynamic_cast<const DTRecHit1D *>(hitp);
377 
378  if (dtrechit) {
379  DTWireId wireid = dtrechit->wireId();
380 
381  if (associatorByWire) {
382  // matching based on DTWireId : take only "valid" SimHits on that wire
383 
384  auto found = mapOfSimHit.find(wireid);
385  if(found != mapOfSimHit.end()) {
386  for (vector<PSimHit_withFlag>::const_iterator hitIT = found->second.begin();
387  hitIT != found->second.end();
388  ++hitIT) {
389 
390  bool valid_hit = hitIT->second;
391  if (valid_hit) simhits.push_back(hitIT->first);
392  }
393  }
394  }
395  else {
396  // matching based on DTDigiSimLink
397 
398  simtrackids = associateDTHitId(dtrechit);
399 
400  for (vector<SimHitIdpr>::const_iterator idIT = simtrackids.begin(); idIT != simtrackids.end(); ++idIT) {
401  uint32_t trId = idIT->first;
402  EncodedEventId evId = idIT->second;
403 
404  auto found = mapOfSimHit.find(wireid);
405  if(found != mapOfSimHit.end()) {
406  for (vector<PSimHit_withFlag>::const_iterator hitIT = found->second.begin();
407  hitIT != found->second.end();
408  ++hitIT) {
409 
410  if (hitIT->first.trackId() == trId &&
411  hitIT->first.eventId() == evId)
412  simhits.push_back(hitIT->first);
413  }
414  }
415  }
416  }
417 
418  } else {
419  edm::LogWarning("DTHitAssociator")<<"*** WARNING in DTHitAssociator::associateHit, null dynamic_cast !";
420  }
421  return simhits;
422 }
423 
425  const PSimHit & simhit) {
426  bool result = true;
427 
428  DTWireId wireid(simhit.detUnitId());
429  const DTLayer* dtlayer = muongeom->layer(wireid.layerId());
430  LocalPoint entryP = simhit.entryPoint();
431  LocalPoint exitP = simhit.exitPoint();
432  const DTTopology &topo = dtlayer->specificTopology();
433  float xwire = topo.wirePosition(wireid.wire());
434  float xEntry = entryP.x() - xwire;
435  float xExit = exitP.x() - xwire;
436  DTTopology::Side entrySide = topo.onWhichBorder(xEntry,entryP.y(),entryP.z());
437  DTTopology::Side exitSide = topo.onWhichBorder(xExit,exitP.y(),exitP.z());
438 
439  bool noParametrisation =
440  (( entrySide == DTTopology::none || exitSide == DTTopology::none ) ||
441  (entrySide == exitSide) ||
442  ((entrySide == DTTopology::xMin && exitSide == DTTopology::xMax) ||
443  (entrySide == DTTopology::xMax && exitSide == DTTopology::xMin)) );
444 
445  // discard hits where parametrization can not be used
446  if (noParametrisation)
447  {
448  result = false;
449  return result;
450  }
451 
452  float x;
453  LocalPoint hitPos = simhit.localPosition();
454 
455  if(fabs(hitPos.z()) < 0.002) {
456  // hit center within 20 um from z=0, no need to extrapolate.
457  x = hitPos.x() - xwire;
458  } else {
459  x = xEntry - (entryP.z()*(xExit-xEntry))/(exitP.z()-entryP.z());
460  }
461 
462  // discard hits where x is out of range of the parametrization (|x|>2.1 cm)
463  x *= 10.; //cm -> mm
464  if (fabs(x) > 21.)
465  result = false;
466 
467  return result;
468 }
float wirePosition(int wireNumber) const
Returns the x position in the layer of a given wire number.
Definition: DTTopology.cc:86
std::vector< SimHitIdpr > associateHitId(const TrackingRecHit &hit) const
std::pair< uint32_t, EncodedEventId > SimHitIdpr
int event() const
get the contents of the subdetector field (should be protected?)
Side onWhichBorder(float x, float y, float z) const
Definition: DTTopology.cc:109
T y() const
Definition: PV3DBase.h:63
RecHitMap mapOfRecHit
edm::InputTag DTdigiTag
std::vector< PSimHit > associateHit(const TrackingRecHit &hit) const
edm::InputTag DTsimhitsTag
int iEvent
Definition: GenABIO.cc:230
Local3DPoint exitPoint() const
Exit point in the local Det frame.
Definition: PSimHit.h:38
edm::InputTag DTrechitTag
SimHitMap mapOfSimHit
Local3DPoint localPosition() const
Definition: PSimHit.h:44
int bunchCrossing() const
get the detector field from this detid
T z() const
Definition: PV3DBase.h:64
tuple result
Definition: query.py:137
float digiTime() const
Return the time (ns) of the digi used to build the rechit.
Definition: DTRecHit1D.h:113
EncodedEventId eventId() const
Definition: PSimHit.h:105
Side
Sides of the cell.
Definition: DTTopology.h:88
bool isValid() const
Definition: HandleBase.h:75
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:413
#define LogTrace(id)
tuple conf
Definition: dbtoconf.py:185
T const * product() const
Definition: Handle.h:81
edm::InputTag DTdigisimlinkTag
const T & get() const
Definition: EventSetup.h:55
std::vector< DTDigi >::const_iterator const_iterator
unsigned short processType() const
Definition: PSimHit.h:118
std::string const & label() const
Definition: InputTag.h:43
#define begin
Definition: vmac.h:30
int particleType() const
Definition: PSimHit.h:85
void initEvent(const edm::Event &, const edm::EventSetup &)
DTHitAssociator(const edm::Event &, const edm::EventSetup &, const edm::ParameterSet &, bool printRtS)
unsigned int trackId() const
Definition: PSimHit.h:102
std::pair< const_iterator, const_iterator > Range
std::vector< PSimHit > PSimHitContainer
edm::InputTag DTsimhitsXFTag
bool SimHitOK(const edm::ESHandle< DTGeometry > &, const PSimHit &)
T x() const
Definition: PV3DBase.h:62
Local3DPoint entryPoint() const
Entry point in the local Det frame.
Definition: PSimHit.h:35
std::vector< SimHitIdpr > associateDTHitId(const DTRecHit1D *dtrechit) const
unsigned int detUnitId() const
Definition: PSimHit.h:93
DTWireId wireId() const
Return the wireId.
Definition: DTRecHit1D.h:107