CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
GlobalRecHitsProducer.cc
Go to the documentation of this file.
1 
14 
16  fName(""), verbosity(0), frequency(0), label(""), getAllProvenances(false),
17  printProvenanceInfo(false), count(0)
18 {
19  std::string MsgLoggerCat = "GlobalRecHitsProducer_GlobalRecHitsProducer";
20 
21  // get information from parameter set
22  fName = iPSet.getUntrackedParameter<std::string>("Name");
23  verbosity = iPSet.getUntrackedParameter<int>("Verbosity");
24  frequency = iPSet.getUntrackedParameter<int>("Frequency");
25  label = iPSet.getParameter<std::string>("Label");
26  edm::ParameterSet m_Prov =
27  iPSet.getParameter<edm::ParameterSet>("ProvenanceLookup");
29  m_Prov.getUntrackedParameter<bool>("GetAllProvenances");
31  m_Prov.getUntrackedParameter<bool>("PrintProvenanceInfo");
32 
33  //get Labels to use to extract information
34  ECalEBSrc_ = iPSet.getParameter<edm::InputTag>("ECalEBSrc");
35  ECalUncalEBSrc_ = iPSet.getParameter<edm::InputTag>("ECalUncalEBSrc");
36  ECalEESrc_ = iPSet.getParameter<edm::InputTag>("ECalEESrc");
37  ECalUncalEESrc_ = iPSet.getParameter<edm::InputTag>("ECalUncalEESrc");
38  ECalESSrc_ = iPSet.getParameter<edm::InputTag>("ECalESSrc");
39  HCalSrc_ = iPSet.getParameter<edm::InputTag>("HCalSrc");
40  SiStripSrc_ = iPSet.getParameter<edm::InputTag>("SiStripSrc");
41  SiPxlSrc_ = iPSet.getParameter<edm::InputTag>("SiPxlSrc");
42  MuDTSrc_ = iPSet.getParameter<edm::InputTag>("MuDTSrc");
43  MuDTSimSrc_ = iPSet.getParameter<edm::InputTag>("MuDTSimSrc");
44  MuCSCSrc_ = iPSet.getParameter<edm::InputTag>("MuCSCSrc");
45  MuRPCSrc_ = iPSet.getParameter<edm::InputTag>("MuRPCSrc");
46  MuRPCSimSrc_ = iPSet.getParameter<edm::InputTag>("MuRPCSimSrc");
47 
48  conf_ = iPSet;
49 
50  // use value of first digit to determine default output level (inclusive)
51  // 0 is none, 1 is basic, 2 is fill output, 3 is gather output
52  verbosity %= 10;
53 
54  // create persistent object
55  produces<PGlobalRecHit>(label);
56 
57  // print out Parameter Set information being used
58  if (verbosity >= 0) {
59  edm::LogInfo(MsgLoggerCat)
60  << "\n===============================\n"
61  << "Initialized as EDProducer with parameter values:\n"
62  << " Name = " << fName << "\n"
63  << " Verbosity = " << verbosity << "\n"
64  << " Frequency = " << frequency << "\n"
65  << " Label = " << label << "\n"
66  << " GetProv = " << getAllProvenances << "\n"
67  << " PrintProv = " << printProvenanceInfo << "\n"
68  << " ECalEBSrc = " << ECalEBSrc_.label()
69  << ":" << ECalEBSrc_.instance() << "\n"
70  << " ECalUncalEBSrc = " << ECalUncalEBSrc_.label()
71  << ":" << ECalUncalEBSrc_.instance() << "\n"
72  << " ECalEESrc = " << ECalEESrc_.label()
73  << ":" << ECalUncalEESrc_.instance() << "\n"
74  << " ECalUncalEESrc = " << ECalUncalEESrc_.label()
75  << ":" << ECalEESrc_.instance() << "\n"
76  << " ECalESSrc = " << ECalESSrc_.label()
77  << ":" << ECalESSrc_.instance() << "\n"
78  << " HCalSrc = " << HCalSrc_.label()
79  << ":" << HCalSrc_.instance() << "\n"
80  << " SiStripSrc = " << SiStripSrc_.label()
81  << ":" << SiStripSrc_.instance() << "\n"
82  << " SiPixelSrc = " << SiPxlSrc_.label()
83  << ":" << SiPxlSrc_.instance() << "\n"
84  << " MuDTSrc = " << MuDTSrc_.label()
85  << ":" << MuDTSrc_.instance() << "\n"
86  << " MuDTSimSrc = " << MuDTSimSrc_.label()
87  << ":" << MuDTSimSrc_.instance() << "\n"
88  << " MuCSCSrc = " << MuCSCSrc_.label()
89  << ":" << MuCSCSrc_.instance() << "\n"
90  << " MuRPCSrc = " << MuRPCSrc_.label()
91  << ":" << MuRPCSrc_.instance() << "\n"
92  << " MuRPCSimSrc = " << MuRPCSimSrc_.label()
93  << ":" << MuRPCSimSrc_.instance() << "\n"
94  << "===============================\n";
95  }
96 }
97 
99 {
100 }
101 
103 {
104  std::string MsgLoggerCat = "GlobalRecHitsProducer_beginJob";
105 
106  // clear storage vectors
107  clear();
108  return;
109 }
110 
112 {
113  std::string MsgLoggerCat = "GlobalRecHitsProducer_endJob";
114  if (verbosity >= 0)
115  edm::LogInfo(MsgLoggerCat)
116  << "Terminating having processed " << count << " events.";
117  return;
118 }
119 
121  const edm::EventSetup& iSetup)
122 {
123  std::string MsgLoggerCat = "GlobalRecHitsProducer_produce";
124 
125  // keep track of number of events processed
126  ++count;
127 
128  // get event id information
129  int nrun = iEvent.id().run();
130  int nevt = iEvent.id().event();
131 
132  if (verbosity > 0) {
133  edm::LogInfo(MsgLoggerCat)
134  << "Processing run " << nrun << ", event " << nevt
135  << " (" << count << " events total)";
136  } else if (verbosity == 0) {
137  if (nevt%frequency == 0 || nevt == 1) {
138  edm::LogInfo(MsgLoggerCat)
139  << "Processing run " << nrun << ", event " << nevt
140  << " (" << count << " events total)";
141  }
142  }
143 
144  // clear event holders
145  clear();
146 
147  // look at information available in the event
148  if (getAllProvenances) {
149 
150  std::vector<const edm::Provenance*> AllProv;
151  iEvent.getAllProvenance(AllProv);
152 
153  if (verbosity >= 0)
154  edm::LogInfo(MsgLoggerCat)
155  << "Number of Provenances = " << AllProv.size();
156 
157  if (printProvenanceInfo && (verbosity >= 0)) {
158  TString eventout("\nProvenance info:\n");
159 
160  for (unsigned int i = 0; i < AllProv.size(); ++i) {
161  eventout += "\n ******************************";
162  eventout += "\n Module : ";
163  //eventout += (AllProv[i]->product).moduleLabel();
164  eventout += AllProv[i]->moduleLabel();
165  eventout += "\n ProductID : ";
166  //eventout += (AllProv[i]->product).productID_.id_;
167  eventout += AllProv[i]->productID().id();
168  eventout += "\n ClassName : ";
169  //eventout += (AllProv[i]->product).fullClassName_;
170  eventout += AllProv[i]->className();
171  eventout += "\n InstanceName : ";
172  //eventout += (AllProv[i]->product).productInstanceName_;
173  eventout += AllProv[i]->productInstanceName();
174  eventout += "\n BranchName : ";
175  //eventout += (AllProv[i]->product).branchName_;
176  eventout += AllProv[i]->branchName();
177  }
178  eventout += "\n ******************************\n";
179  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
180  printProvenanceInfo = false;
181  }
182  getAllProvenances = false;
183  }
184 
185  // call fill functions
186  // gather Ecal information from event
187  fillECal(iEvent, iSetup);
188  // gather Hcal information from event
189  fillHCal(iEvent, iSetup);
190  // gather Track information from event
191  fillTrk(iEvent, iSetup);
192  // gather Muon information from event
193  fillMuon(iEvent, iSetup);
194 
195  if (verbosity > 0)
196  edm::LogInfo (MsgLoggerCat)
197  << "Done gathering data from event.";
198 
199  // produce object to put into event
200  std::auto_ptr<PGlobalRecHit> pOut(new PGlobalRecHit);
201 
202  if (verbosity > 2)
203  edm::LogInfo (MsgLoggerCat)
204  << "Saving event contents:";
205 
206  // call store functions
207  // store ECal information in produce
208  storeECal(*pOut);
209  // store HCal information in produce
210  storeHCal(*pOut);
211  // store Track information in produce
212  storeTrk(*pOut);
213  // store Muon information in produce
214  storeMuon(*pOut);
215 
216  // store information in event
217  iEvent.put(pOut,label);
218 
219  return;
220 }
221 
223  const edm::EventSetup& iSetup)
224 {
225  std::string MsgLoggerCat = "GlobalRecHitsProducer_fillECal";
226 
227  TString eventout;
228  if (verbosity > 0)
229  eventout = "\nGathering info:";
230 
231  // extract crossing frame from event
232  //edm::Handle<CrossingFrame> crossingFrame;
233  edm::Handle<CrossingFrame<PCaloHit> > crossingFrame;
234  //iEvent.getByType(crossingFrame);
235  //if (!crossingFrame.isValid()) {
236  // edm::LogWarning(MsgLoggerCat)
237  // << "Unable to find crossingFrame in event!";
238  // return;
239  //}
240 
242  //extract EB information
245  iEvent.getByLabel(ECalUncalEBSrc_, EcalUncalibRecHitEB);
246  if (!EcalUncalibRecHitEB.isValid()) {
247  edm::LogWarning(MsgLoggerCat)
248  << "Unable to find EcalUncalRecHitEB in event!";
249  return;
250  }
251 
252  edm::Handle<EBRecHitCollection> EcalRecHitEB;
253  iEvent.getByLabel(ECalEBSrc_, EcalRecHitEB);
254  if (!EcalRecHitEB.isValid()) {
255  edm::LogWarning(MsgLoggerCat)
256  << "Unable to find EcalRecHitEB in event!";
257  return;
258  }
259 
260  // loop over simhits
261  const std::string barrelHitsName("EcalHitsEB");
262  iEvent.getByLabel("mix",barrelHitsName,crossingFrame);
263  if (!crossingFrame.isValid()) {
264  edm::LogWarning(MsgLoggerCat)
265  << "Unable to find cal barrel crossingFrame in event!";
266  return;
267  }
268  //std::auto_ptr<MixCollection<PCaloHit> >
269  // barrelHits(new MixCollection<PCaloHit>
270  // (crossingFrame.product(), barrelHitsName));
271  std::auto_ptr<MixCollection<PCaloHit> >
272  barrelHits(new MixCollection<PCaloHit>(crossingFrame.product()));
273 
274  // keep track of sum of simhit energy in each crystal
275  MapType ebSimMap;
277  = barrelHits->begin();
278  hitItr != barrelHits->end();
279  ++hitItr) {
280 
281  EBDetId ebid = EBDetId(hitItr->id());
282 
283  uint32_t crystid = ebid.rawId();
284  ebSimMap[crystid] += hitItr->energy();
285  }
286 
287  int nEBRecHits = 0;
288  // loop over RecHits
289  const EBUncalibratedRecHitCollection *EBUncalibRecHit =
290  EcalUncalibRecHitEB.product();
291  const EBRecHitCollection *EBRecHit = EcalRecHitEB.product();
292 
294  EBUncalibRecHit->begin();
295  uncalibRecHit != EBUncalibRecHit->end();
296  ++uncalibRecHit) {
297 
298  EBDetId EBid = EBDetId(uncalibRecHit->id());
299 
300  EcalRecHitCollection::const_iterator myRecHit = EBRecHit->find(EBid);
301 
302  if (myRecHit != EBRecHit->end()) {
303  ++nEBRecHits;
304  EBRE.push_back(myRecHit->energy());
305  EBSHE.push_back(ebSimMap[EBid.rawId()]);
306  }
307  }
308 
309  if (verbosity > 1) {
310  eventout += "\n Number of EBRecHits collected:............ ";
311  eventout += nEBRecHits;
312  }
313 
315  //extract EE information
318  iEvent.getByLabel(ECalUncalEESrc_, EcalUncalibRecHitEE);
319  if (!EcalUncalibRecHitEE.isValid()) {
320  edm::LogWarning(MsgLoggerCat)
321  << "Unable to find EcalUncalRecHitEE in event!";
322  return;
323  }
324 
325  edm::Handle<EERecHitCollection> EcalRecHitEE;
326  iEvent.getByLabel(ECalEESrc_, EcalRecHitEE);
327  if (!EcalRecHitEE.isValid()) {
328  edm::LogWarning(MsgLoggerCat)
329  << "Unable to find EcalRecHitEE in event!";
330  return;
331  }
332 
333  // loop over simhits
334  const std::string endcapHitsName("EcalHitsEE");
335  iEvent.getByLabel("mix",endcapHitsName,crossingFrame);
336  if (!crossingFrame.isValid()) {
337  edm::LogWarning(MsgLoggerCat)
338  << "Unable to find cal endcap crossingFrame in event!";
339  return;
340  }
341  //std::auto_ptr<MixCollection<PCaloHit> >
342  // endcapHits(new MixCollection<PCaloHit>
343  // (crossingFrame.product(), endcapHitsName));
344  std::auto_ptr<MixCollection<PCaloHit> >
345  endcapHits(new MixCollection<PCaloHit>(crossingFrame.product()));
346 
347  // keep track of sum of simhit energy in each crystal
348  MapType eeSimMap;
350  = endcapHits->begin();
351  hitItr != endcapHits->end();
352  ++hitItr) {
353 
354  EEDetId eeid = EEDetId(hitItr->id());
355 
356  uint32_t crystid = eeid.rawId();
357  eeSimMap[crystid] += hitItr->energy();
358  }
359 
360  int nEERecHits = 0;
361  // loop over RecHits
362  const EEUncalibratedRecHitCollection *EEUncalibRecHit =
363  EcalUncalibRecHitEE.product();
364  const EERecHitCollection *EERecHit = EcalRecHitEE.product();
365 
367  EEUncalibRecHit->begin();
368  uncalibRecHit != EEUncalibRecHit->end();
369  ++uncalibRecHit) {
370 
371  EEDetId EEid = EEDetId(uncalibRecHit->id());
372 
373  EcalRecHitCollection::const_iterator myRecHit = EERecHit->find(EEid);
374 
375  if (myRecHit != EERecHit->end()) {
376  ++nEERecHits;
377  EERE.push_back(myRecHit->energy());
378  EESHE.push_back(eeSimMap[EEid.rawId()]);
379  }
380  }
381 
382  if (verbosity > 1) {
383  eventout += "\n Number of EERecHits collected:............ ";
384  eventout += nEERecHits;
385  }
386 
388  //extract ES information
390  edm::Handle<ESRecHitCollection> EcalRecHitES;
391  iEvent.getByLabel(ECalESSrc_, EcalRecHitES);
392  if (!EcalRecHitES.isValid()) {
393  edm::LogWarning(MsgLoggerCat)
394  << "Unable to find EcalRecHitES in event!";
395  return;
396  }
397 
398  // loop over simhits
399  const std::string preshowerHitsName("EcalHitsES");
400  iEvent.getByLabel("mix",preshowerHitsName,crossingFrame);
401  if (!crossingFrame.isValid()) {
402  edm::LogWarning(MsgLoggerCat)
403  << "Unable to find cal preshower crossingFrame in event!";
404  return;
405  }
406  //std::auto_ptr<MixCollection<PCaloHit> >
407  // preshowerHits(new MixCollection<PCaloHit>
408  // (crossingFrame.product(), preshowerHitsName));
409  std::auto_ptr<MixCollection<PCaloHit> >
410  preshowerHits(new MixCollection<PCaloHit>(crossingFrame.product()));
411 
412  // keep track of sum of simhit energy in each crystal
413  MapType esSimMap;
415  = preshowerHits->begin();
416  hitItr != preshowerHits->end();
417  ++hitItr) {
418 
419  ESDetId esid = ESDetId(hitItr->id());
420 
421  uint32_t crystid = esid.rawId();
422  esSimMap[crystid] += hitItr->energy();
423  }
424 
425  int nESRecHits = 0;
426  // loop over RecHits
427  const ESRecHitCollection *ESRecHit = EcalRecHitES.product();
429  ESRecHit->begin();
430  recHit != ESRecHit->end();
431  ++recHit) {
432 
433  ESDetId ESid = ESDetId(recHit->id());
434 
435  ++nESRecHits;
436  ESRE.push_back(recHit->energy());
437  ESSHE.push_back(esSimMap[ESid.rawId()]);
438  }
439 
440  if (verbosity > 1) {
441  eventout += "\n Number of ESRecHits collected:............ ";
442  eventout += nESRecHits;
443  }
444 
445  if (verbosity > 0)
446  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
447 
448  return;
449 }
450 
452 {
453  std::string MsgLoggerCat = "GlobalRecHitsProducer_storeECal";
454 
455  if (verbosity > 2) {
456  TString eventout("\n nEBRecHits = ");
457  eventout += EBRE.size();
458  for (unsigned int i = 0; i < EBRE.size(); ++i) {
459  eventout += "\n (RE, SHE) = (";
460  eventout += EBRE[i];
461  eventout += ", ";
462  eventout += EBSHE[i];
463  eventout += ")";
464  }
465  eventout += "\n nEERecHits = ";
466  eventout += EERE.size();
467  for (unsigned int i = 0; i < EERE.size(); ++i) {
468  eventout += "\n (RE, SHE) = (";
469  eventout += EERE[i];
470  eventout += ", ";
471  eventout += EESHE[i];
472  eventout += ")";
473  }
474  eventout += "\n nESRecHits = ";
475  eventout += ESRE.size();
476  for (unsigned int i = 0; i < ESRE.size(); ++i) {
477  eventout += "\n (RE, SHE) = (";
478  eventout += ESRE[i];
479  eventout += ", ";
480  eventout += ESSHE[i];
481  eventout += ")";
482  }
483  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
484  }
485 
486  product.putEBCalRecHits(EBRE,EBSHE);
487  product.putEECalRecHits(EERE,EESHE);
488  product.putESCalRecHits(ESRE,ESSHE);
489 
490  return;
491 }
492 
494  const edm::EventSetup& iSetup)
495 {
496  std::string MsgLoggerCat = "GlobalRecHitsProducer_fillHCal";
497 
498  TString eventout;
499  if (verbosity > 0)
500  eventout = "\nGathering info:";
501 
502  // get geometry
504  iSetup.get<CaloGeometryRecord>().get(geometry);
505  if (!geometry.isValid()) {
506  edm::LogWarning(MsgLoggerCat)
507  << "Unable to find CaloGeometry in event!";
508  return;
509  }
510 
512  // extract simhit info
515  iEvent.getByLabel(HCalSrc_,hcalHits);
516  if (!hcalHits.isValid()) {
517  edm::LogWarning(MsgLoggerCat)
518  << "Unable to find hcalHits in event!";
519  return;
520  }
521  const edm::PCaloHitContainer *simhitResult = hcalHits.product();
522 
523  MapType fHBEnergySimHits;
524  MapType fHEEnergySimHits;
525  MapType fHOEnergySimHits;
526  MapType fHFEnergySimHits;
527  for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin();
528  simhits != simhitResult->end();
529  ++simhits) {
530 
531  HcalDetId detId(simhits->id());
532  uint32_t cellid = detId.rawId();
533 
534  if (detId.subdet() == sdHcalBrl){
535  fHBEnergySimHits[cellid] += simhits->energy();
536  }
537  if (detId.subdet() == sdHcalEC){
538  fHEEnergySimHits[cellid] += simhits->energy();
539  }
540  if (detId.subdet() == sdHcalOut){
541  fHOEnergySimHits[cellid] += simhits->energy();
542  }
543  if (detId.subdet() == sdHcalFwd){
544  fHFEnergySimHits[cellid] += simhits->energy();
545  }
546  }
547 
548  // max values to be used (HO is found in HB)
549  Double_t maxHBEnergy = 0.;
550  Double_t maxHEEnergy = 0.;
551  Double_t maxHFEnergy = 0.;
552 
553  Double_t maxHBPhi = -1000.;
554  Double_t maxHEPhi = -1000.;
555  Double_t maxHOPhi = -1000.;
556  Double_t maxHFPhi = -1000.;
557 
558  Double_t maxHBEta = -1000.;
559  Double_t maxHEEta = -1000.;
560  Double_t maxHOEta = -1000.;
561  Double_t maxHFEta = -1000.;
562 
563  Double_t PI = 3.141592653589;
564 
566  // get HBHE information
568  std::vector<edm::Handle<HBHERecHitCollection> > hbhe;
569  iEvent.getManyByType(hbhe);
570  if (!hbhe[0].isValid()) {
571  edm::LogWarning(MsgLoggerCat)
572  << "Unable to find any HBHERecHitCollections in event!";
573  return;
574  }
575  std::vector<edm::Handle<HBHERecHitCollection> >::iterator ihbhe;
576 
577  int iHB = 0;
578  int iHE = 0;
579  for (ihbhe = hbhe.begin(); ihbhe != hbhe.end(); ++ihbhe) {
580 
581  // find max values
582  for (HBHERecHitCollection::const_iterator jhbhe = (*ihbhe)->begin();
583  jhbhe != (*ihbhe)->end(); ++jhbhe) {
584 
585  HcalDetId cell(jhbhe->id());
586 
587  if (cell.subdet() == sdHcalBrl) {
588 
589  const CaloCellGeometry* cellGeometry =
590  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
591  double fEta = cellGeometry->getPosition().eta () ;
592  double fPhi = cellGeometry->getPosition().phi () ;
593  if ( (jhbhe->energy()) > maxHBEnergy ) {
594  maxHBEnergy = jhbhe->energy();
595  maxHBPhi = fPhi;
596  maxHOPhi = maxHBPhi;
597  maxHBEta = fEta;
598  maxHOEta = maxHBEta;
599  }
600  }
601 
602  if (cell.subdet() == sdHcalEC) {
603 
604  const CaloCellGeometry* cellGeometry =
605  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
606  double fEta = cellGeometry->getPosition().eta () ;
607  double fPhi = cellGeometry->getPosition().phi () ;
608  if ( (jhbhe->energy()) > maxHEEnergy ) {
609  maxHEEnergy = jhbhe->energy();
610  maxHEPhi = fPhi;
611  maxHEEta = fEta;
612  }
613  }
614  } // end find max values
615 
616  for (HBHERecHitCollection::const_iterator jhbhe = (*ihbhe)->begin();
617  jhbhe != (*ihbhe)->end(); ++jhbhe) {
618 
619  HcalDetId cell(jhbhe->id());
620 
621  if (cell.subdet() == sdHcalBrl) {
622 
623  ++iHB;
624 
625  const CaloCellGeometry* cellGeometry =
626  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
627  double fEta = cellGeometry->getPosition().eta () ;
628  double fPhi = cellGeometry->getPosition().phi () ;
629 
630  float deltaphi = maxHBPhi - fPhi;
631  if (fPhi > maxHBPhi) { deltaphi = fPhi - maxHBPhi;}
632  if (deltaphi > PI) { deltaphi = 2.0 * PI - deltaphi;}
633  float deltaeta = fEta - maxHBEta;
634  Double_t r = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
635 
636  HBCalREC.push_back(jhbhe->energy());
637  HBCalR.push_back(r);
638  HBCalSHE.push_back(fHBEnergySimHits[cell.rawId()]);
639  }
640 
641  if (cell.subdet() == sdHcalEC) {
642 
643  ++iHE;
644 
645  const CaloCellGeometry* cellGeometry =
646  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
647  double fEta = cellGeometry->getPosition().eta () ;
648  double fPhi = cellGeometry->getPosition().phi () ;
649 
650  float deltaphi = maxHEPhi - fPhi;
651  if (fPhi > maxHEPhi) { deltaphi = fPhi - maxHEPhi;}
652  if (deltaphi > PI) { deltaphi = 2.0 * PI - deltaphi;}
653  float deltaeta = fEta - maxHEEta;
654  Double_t r = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
655 
656  HECalREC.push_back(jhbhe->energy());
657  HECalR.push_back(r);
658  HECalSHE.push_back(fHEEnergySimHits[cell.rawId()]);
659  }
660  }
661  } // end loop through collection
662 
663 
664  if (verbosity > 1) {
665  eventout += "\n Number of HBRecHits collected:............ ";
666  eventout += iHB;
667  }
668 
669  if (verbosity > 1) {
670  eventout += "\n Number of HERecHits collected:............ ";
671  eventout += iHE;
672  }
673 
675  // get HF information
677  std::vector<edm::Handle<HFRecHitCollection> > hf;
678  iEvent.getManyByType(hf);
679  if (!hf[0].isValid()) {
680  edm::LogWarning(MsgLoggerCat)
681  << "Unable to find any HFRecHitCollections in event!";
682  return;
683  }
684  std::vector<edm::Handle<HFRecHitCollection> >::iterator ihf;
685 
686  int iHF = 0;
687  for (ihf = hf.begin(); ihf != hf.end(); ++ihf) {
688 
689  // find max values
690  for (HFRecHitCollection::const_iterator jhf = (*ihf)->begin();
691  jhf != (*ihf)->end(); ++jhf) {
692 
693  HcalDetId cell(jhf->id());
694 
695  if (cell.subdet() == sdHcalFwd) {
696 
697  const CaloCellGeometry* cellGeometry =
698  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
699  double fEta = cellGeometry->getPosition().eta () ;
700  double fPhi = cellGeometry->getPosition().phi () ;
701  if ( (jhf->energy()) > maxHFEnergy ) {
702  maxHFEnergy = jhf->energy();
703  maxHFPhi = fPhi;
704  maxHFEta = fEta;
705  }
706  }
707  } // end find max values
708 
709  for (HFRecHitCollection::const_iterator jhf = (*ihf)->begin();
710  jhf != (*ihf)->end(); ++jhf) {
711 
712  HcalDetId cell(jhf->id());
713 
714  if (cell.subdet() == sdHcalFwd) {
715 
716  ++iHF;
717 
718  const CaloCellGeometry* cellGeometry =
719  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
720  double fEta = cellGeometry->getPosition().eta () ;
721  double fPhi = cellGeometry->getPosition().phi () ;
722 
723  float deltaphi = maxHBPhi - fPhi;
724  if (fPhi > maxHFPhi) { deltaphi = fPhi - maxHFPhi;}
725  if (deltaphi > PI) { deltaphi = 2.0 * PI - deltaphi;}
726  float deltaeta = fEta - maxHFEta;
727  Double_t r = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
728 
729  HFCalREC.push_back(jhf->energy());
730  HFCalR.push_back(r);
731  HFCalSHE.push_back(fHFEnergySimHits[cell.rawId()]);
732  }
733  }
734  } // end loop through collection
735 
736  if (verbosity > 1) {
737  eventout += "\n Number of HFDigis collected:.............. ";
738  eventout += iHF;
739  }
740 
742  // get HO information
744  std::vector<edm::Handle<HORecHitCollection> > ho;
745  iEvent.getManyByType(ho);
746  if (!ho[0].isValid()) {
747  edm::LogWarning(MsgLoggerCat)
748  << "Unable to find any HORecHitCollections in event!";
749  return;
750  }
751  std::vector<edm::Handle<HORecHitCollection> >::iterator iho;
752 
753  int iHO = 0;
754  for (iho = ho.begin(); iho != ho.end(); ++iho) {
755 
756  for (HORecHitCollection::const_iterator jho = (*iho)->begin();
757  jho != (*iho)->end(); ++jho) {
758 
759  HcalDetId cell(jho->id());
760 
761  if (cell.subdet() == sdHcalOut) {
762 
763  ++iHO;
764 
765  const CaloCellGeometry* cellGeometry =
766  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
767  double fEta = cellGeometry->getPosition().eta () ;
768  double fPhi = cellGeometry->getPosition().phi () ;
769 
770  float deltaphi = maxHOPhi - fPhi;
771  if (fPhi > maxHOPhi) { deltaphi = fPhi - maxHOPhi;}
772  if (deltaphi > PI) { deltaphi = 2.0 * PI - deltaphi;}
773  float deltaeta = fEta - maxHOEta;
774  Double_t r = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
775 
776  HOCalREC.push_back(jho->energy());
777  HOCalR.push_back(r);
778  HOCalSHE.push_back(fHOEnergySimHits[cell.rawId()]);
779  }
780  }
781  } // end loop through collection
782 
783  if (verbosity > 1) {
784  eventout += "\n Number of HODigis collected:.............. ";
785  eventout += iHO;
786  }
787 
788  if (verbosity > 0)
789  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
790 
791  return;
792 }
793 
795 {
796  std::string MsgLoggerCat = "GlobalRecHitsProducer_storeHCal";
797 
798  if (verbosity > 2) {
799  TString eventout("\n nHBRecHits = ");
800  eventout += HBCalREC.size();
801  for (unsigned int i = 0; i < HBCalREC.size(); ++i) {
802  eventout += "\n (REC, R, SHE) = (";
803  eventout += HBCalREC[i];
804  eventout += ", ";
805  eventout += HBCalR[i];
806  eventout += ", ";
807  eventout += HBCalSHE[i];
808  eventout += ")";
809  }
810  eventout += "\n nHERecHits = ";
811  eventout += HECalREC.size();
812  for (unsigned int i = 0; i < HECalREC.size(); ++i) {
813  eventout += "\n (REC, R, SHE) = (";
814  eventout += HECalREC[i];
815  eventout += ", ";
816  eventout += HECalR[i];
817  eventout += ", ";
818  eventout += HECalSHE[i];
819  eventout += ")";
820  }
821  eventout += "\n nHFRecHits = ";
822  eventout += HFCalREC.size();
823  for (unsigned int i = 0; i < HFCalREC.size(); ++i) {
824  eventout += "\n (REC, R, SHE) = (";
825  eventout += HFCalREC[i];
826  eventout += ", ";
827  eventout += HFCalR[i];
828  eventout += ", ";
829  eventout += HFCalSHE[i];
830  eventout += ")";
831  }
832  eventout += "\n nHORecHits = ";
833  eventout += HOCalREC.size();
834  for (unsigned int i = 0; i < HOCalREC.size(); ++i) {
835  eventout += "\n (REC, R, SHE) = (";
836  eventout += HOCalREC[i];
837  eventout += ", ";
838  eventout += HOCalR[i];
839  eventout += ", ";
840  eventout += HOCalSHE[i];
841  eventout += ")";
842  }
843 
844  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
845  }
846 
851 
852  return;
853 }
854 
856  const edm::EventSetup& iSetup)
857 {
858  //Retrieve tracker topology from geometry
859  edm::ESHandle<TrackerTopology> tTopoHandle;
860  iSetup.get<IdealGeometryRecord>().get(tTopoHandle);
861  const TrackerTopology* const tTopo = tTopoHandle.product();
862 
863 
864  std::string MsgLoggerCat = "GlobalRecHitsProducer_fillTrk";
865 
866  TString eventout;
867  if (verbosity > 0)
868  eventout = "\nGathering info:";
869 
870  // get strip information
872  iEvent.getByLabel(SiStripSrc_, rechitsmatched);
873  if (!rechitsmatched.isValid()) {
874  edm::LogWarning(MsgLoggerCat)
875  << "Unable to find stripmatchedrechits in event!";
876  return;
877  }
878 
879  TrackerHitAssociator associate(iEvent,conf_);
880 
882  iSetup.get<TrackerDigiGeometryRecord>().get(pDD);
883  if (!pDD.isValid()) {
884  edm::LogWarning(MsgLoggerCat)
885  << "Unable to find TrackerDigiGeometry in event!";
886  return;
887  }
888  const TrackerGeometry &tracker(*pDD);
889 
890  int nStripBrl = 0, nStripFwd = 0;
891 
892  // loop over det units
893  for (TrackerGeometry::DetContainer::const_iterator it = pDD->dets().begin();
894  it != pDD->dets().end(); ++it) {
895 
896  uint32_t myid = ((*it)->geographicalId()).rawId();
897  DetId detid = ((*it)->geographicalId());
898 
899  //loop over rechits-matched in the same subdetector
900  SiStripMatchedRecHit2DCollection::const_iterator rechitmatchedMatch = rechitsmatched->find(detid);
901 
902  if (rechitmatchedMatch != rechitsmatched->end()) {
903  SiStripMatchedRecHit2DCollection::DetSet rechitmatchedRange = *rechitmatchedMatch;
904  SiStripMatchedRecHit2DCollection::DetSet::const_iterator rechitmatchedRangeIteratorBegin = rechitmatchedRange.begin();
905  SiStripMatchedRecHit2DCollection::DetSet::const_iterator rechitmatchedRangeIteratorEnd = rechitmatchedRange.end();
906  SiStripMatchedRecHit2DCollection::DetSet::const_iterator itermatched = rechitmatchedRangeIteratorBegin;
907 
908  for ( itermatched = rechitmatchedRangeIteratorBegin;
909  itermatched != rechitmatchedRangeIteratorEnd;
910  ++itermatched) {
911 
912  SiStripMatchedRecHit2D const rechit = *itermatched;
913  LocalPoint position = rechit.localPosition();
914 
915  float mindist = 999999.;
916  float distx = 999999.;
917  float disty = 999999.;
918  float dist = 999999.;
919  std::pair<LocalPoint,LocalVector> closestPair;
920  matched.clear();
921 
922  float rechitmatchedx = position.x();
923  float rechitmatchedy = position.y();
924 
925  matched = associate.associateHit(rechit);
926 
927  if (!matched.empty()) {
928  //project simhit;
929  const GluedGeomDet* gluedDet =
930  (const GluedGeomDet*)tracker.idToDet(rechit.geographicalId());
931  const StripGeomDetUnit* partnerstripdet =
932  (StripGeomDetUnit*) gluedDet->stereoDet();
933  std::pair<LocalPoint,LocalVector> hitPair;
934 
935  for(std::vector<PSimHit>::const_iterator m = matched.begin();
936  m != matched.end(); m++){
937  //project simhit;
938  hitPair = projectHit((*m),partnerstripdet,gluedDet->surface());
939  distx = fabs(rechitmatchedx - hitPair.first.x());
940  disty = fabs(rechitmatchedy - hitPair.first.y());
941  dist = sqrt(distx*distx+disty*disty);
942 
943  if(dist < mindist){
944  mindist = dist;
945  closestPair = hitPair;
946  }
947  }
948 
949  // get TIB
950  if (detid.subdetId() == sdSiTIB) {
951 
952 
953  ++nStripBrl;
954 
955  if (tTopo->tibLayer(myid) == 1) {
956  TIBL1RX.push_back(rechitmatchedx);
957  TIBL1RY.push_back(rechitmatchedy);
958  TIBL1SX.push_back(closestPair.first.x());
959  TIBL1SY.push_back(closestPair.first.y());
960  }
961  if (tTopo->tibLayer(myid) == 2) {
962  TIBL2RX.push_back(rechitmatchedx);
963  TIBL2RY.push_back(rechitmatchedy);
964  TIBL2SX.push_back(closestPair.first.x());
965  TIBL2SY.push_back(closestPair.first.y());
966  }
967  if (tTopo->tibLayer(myid) == 3) {
968  TIBL3RX.push_back(rechitmatchedx);
969  TIBL3RY.push_back(rechitmatchedy);
970  TIBL3SX.push_back(closestPair.first.x());
971  TIBL3SY.push_back(closestPair.first.y());
972  }
973  if (tTopo->tibLayer(myid) == 4) {
974  TIBL4RX.push_back(rechitmatchedx);
975  TIBL4RY.push_back(rechitmatchedy);
976  TIBL4SX.push_back(closestPair.first.x());
977  TIBL4SY.push_back(closestPair.first.y());
978  }
979  }
980 
981  // get TOB
982  if (detid.subdetId() == sdSiTOB) {
983 
984 
985  ++nStripBrl;
986 
987  if (tTopo->tobLayer(myid) == 1) {
988  TOBL1RX.push_back(rechitmatchedx);
989  TOBL1RY.push_back(rechitmatchedy);
990  TOBL1SX.push_back(closestPair.first.x());
991  TOBL1SY.push_back(closestPair.first.y());
992  }
993  if (tTopo->tobLayer(myid) == 2) {
994  TOBL2RX.push_back(rechitmatchedx);
995  TOBL2RY.push_back(rechitmatchedy);
996  TOBL2SX.push_back(closestPair.first.x());
997  TOBL2SY.push_back(closestPair.first.y());
998  }
999  if (tTopo->tobLayer(myid) == 3) {
1000  TOBL3RX.push_back(rechitmatchedx);
1001  TOBL3RY.push_back(rechitmatchedy);
1002  TOBL3SX.push_back(closestPair.first.x());
1003  TOBL3SY.push_back(closestPair.first.y());
1004  }
1005  if (tTopo->tobLayer(myid) == 4) {
1006  TOBL4RX.push_back(rechitmatchedx);
1007  TOBL4RY.push_back(rechitmatchedy);
1008  TOBL4SX.push_back(closestPair.first.x());
1009  TOBL4SY.push_back(closestPair.first.y());
1010  }
1011  }
1012 
1013  // get TID
1014  if (detid.subdetId() == sdSiTID) {
1015 
1016 
1017  ++nStripFwd;
1018 
1019  if (tTopo->tidWheel(myid) == 1) {
1020  TIDW1RX.push_back(rechitmatchedx);
1021  TIDW1RY.push_back(rechitmatchedy);
1022  TIDW1SX.push_back(closestPair.first.x());
1023  TIDW1SY.push_back(closestPair.first.y());
1024  }
1025  if (tTopo->tidWheel(myid) == 2) {
1026  TIDW2RX.push_back(rechitmatchedx);
1027  TIDW2RY.push_back(rechitmatchedy);
1028  TIDW2SX.push_back(closestPair.first.x());
1029  TIDW2SY.push_back(closestPair.first.y());
1030  }
1031  if (tTopo->tidWheel(myid) == 3) {
1032  TIDW3RX.push_back(rechitmatchedx);
1033  TIDW3RY.push_back(rechitmatchedy);
1034  TIDW3SX.push_back(closestPair.first.x());
1035  TIDW3SY.push_back(closestPair.first.y());
1036  }
1037  }
1038 
1039  // get TEC
1040  if (detid.subdetId() == sdSiTEC) {
1041 
1042 
1043  ++nStripFwd;
1044 
1045  if (tTopo->tecWheel(myid) == 1) {
1046  TECW1RX.push_back(rechitmatchedx);
1047  TECW1RY.push_back(rechitmatchedy);
1048  TECW1SX.push_back(closestPair.first.x());
1049  TECW1SY.push_back(closestPair.first.y());
1050  }
1051  if (tTopo->tecWheel(myid) == 2) {
1052  TECW2RX.push_back(rechitmatchedx);
1053  TECW2RY.push_back(rechitmatchedy);
1054  TECW2SX.push_back(closestPair.first.x());
1055  TECW2SY.push_back(closestPair.first.y());
1056  }
1057  if (tTopo->tecWheel(myid) == 3) {
1058  TECW3RX.push_back(rechitmatchedx);
1059  TECW3RY.push_back(rechitmatchedy);
1060  TECW3SX.push_back(closestPair.first.x());
1061  TECW3SY.push_back(closestPair.first.y());
1062  }
1063  if (tTopo->tecWheel(myid) == 4) {
1064  TECW4RX.push_back(rechitmatchedx);
1065  TECW4RY.push_back(rechitmatchedy);
1066  TECW4SX.push_back(closestPair.first.x());
1067  TECW4SY.push_back(closestPair.first.y());
1068  }
1069  if (tTopo->tecWheel(myid) == 5) {
1070  TECW5RX.push_back(rechitmatchedx);
1071  TECW5RY.push_back(rechitmatchedy);
1072  TECW5SX.push_back(closestPair.first.x());
1073  TECW5SY.push_back(closestPair.first.y());
1074  }
1075  if (tTopo->tecWheel(myid) == 6) {
1076  TECW6RX.push_back(rechitmatchedx);
1077  TECW6RY.push_back(rechitmatchedy);
1078  TECW6SX.push_back(closestPair.first.x());
1079  TECW6SY.push_back(closestPair.first.y());
1080  }
1081  if (tTopo->tecWheel(myid) == 7) {
1082  TECW7RX.push_back(rechitmatchedx);
1083  TECW7RY.push_back(rechitmatchedy);
1084  TECW7SX.push_back(closestPair.first.x());
1085  TECW7SY.push_back(closestPair.first.y());
1086  }
1087  if (tTopo->tecWheel(myid) == 8) {
1088  TECW8RX.push_back(rechitmatchedx);
1089  TECW8RY.push_back(rechitmatchedy);
1090  TECW8SX.push_back(closestPair.first.x());
1091  TECW8SY.push_back(closestPair.first.y());
1092  }
1093  }
1094 
1095  } // end if matched empty
1096  }
1097  }
1098  } // end loop over det units
1099 
1100  if (verbosity > 1) {
1101  eventout += "\n Number of BrlStripRecHits collected:...... ";
1102  eventout += nStripBrl;
1103  }
1104 
1105  if (verbosity > 1) {
1106  eventout += "\n Number of FrwdStripRecHits collected:..... ";
1107  eventout += nStripFwd;
1108  }
1109 
1110  // get pixel information
1111  //Get RecHits
1113  iEvent.getByLabel(SiPxlSrc_, recHitColl);
1114  if (!recHitColl.isValid()) {
1115  edm::LogWarning(MsgLoggerCat)
1116  << "Unable to find SiPixelRecHitCollection in event!";
1117  return;
1118  }
1119 
1120  //Get event setup
1122  iSetup.get<TrackerDigiGeometryRecord>().get(geom);
1123  if (!geom.isValid()) {
1124  edm::LogWarning(MsgLoggerCat)
1125  << "Unable to find TrackerDigiGeometry in event!";
1126  return;
1127  }
1128  //const TrackerGeometry& theTracker(*geom);
1129 
1130  int nPxlBrl = 0, nPxlFwd = 0;
1131  //iterate over detunits
1132  for (TrackerGeometry::DetContainer::const_iterator it = geom->dets().begin();
1133  it != geom->dets().end(); ++it) {
1134 
1135  uint32_t myid = ((*it)->geographicalId()).rawId();
1136  DetId detId = ((*it)->geographicalId());
1137  int subid = detId.subdetId();
1138 
1139  if (! ((subid == sdPxlBrl) || (subid == sdPxlFwd))) continue;
1140 
1141  //const PixelGeomDetUnit * theGeomDet =
1142  // dynamic_cast<const PixelGeomDetUnit*>(theTracker.idToDet(detId) );
1143 
1144  SiPixelRecHitCollection::const_iterator pixeldet = recHitColl->find(detId);
1145  if (pixeldet == recHitColl->end()) continue;
1146  SiPixelRecHitCollection::DetSet pixelrechitRange = *pixeldet;
1147  SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorBegin = pixelrechitRange.begin();
1148  SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorEnd = pixelrechitRange.end();
1149  SiPixelRecHitCollection::DetSet::const_iterator pixeliter = pixelrechitRangeIteratorBegin;
1150  std::vector<PSimHit> matched;
1151 
1152  //----Loop over rechits for this detId
1153  for ( ; pixeliter != pixelrechitRangeIteratorEnd; ++pixeliter) {
1154 
1155  matched.clear();
1156  matched = associate.associateHit(*pixeliter);
1157 
1158  if ( !matched.empty() ) {
1159 
1160  float closest = 9999.9;
1161  //std::vector<PSimHit>::const_iterator closestit = matched.begin();
1162  LocalPoint lp = pixeliter->localPosition();
1163  float rechit_x = lp.x();
1164  float rechit_y = lp.y();
1165 
1166  float sim_x = 0.;
1167  float sim_y = 0.;
1168 
1169  //loop over sim hits and fill closet
1170  for (std::vector<PSimHit>::const_iterator m = matched.begin();
1171  m != matched.end(); ++m) {
1172 
1173  float sim_x1 = (*m).entryPoint().x();
1174  float sim_x2 = (*m).exitPoint().x();
1175  float sim_xpos = 0.5*(sim_x1+sim_x2);
1176 
1177  float sim_y1 = (*m).entryPoint().y();
1178  float sim_y2 = (*m).exitPoint().y();
1179  float sim_ypos = 0.5*(sim_y1+sim_y2);
1180 
1181  float x_res = fabs(sim_xpos - rechit_x);
1182  float y_res = fabs(sim_ypos - rechit_y);
1183 
1184  float dist = sqrt(x_res*x_res + y_res*y_res);
1185 
1186  if ( dist < closest ) {
1187  closest = dist;
1188  sim_x = sim_xpos;
1189  sim_y = sim_ypos;
1190  }
1191  } // end sim hit loop
1192 
1193  // get Barrel pixels
1194  if (subid == sdPxlBrl) {
1195 
1196  ++nPxlBrl;
1197 
1198  if (tTopo->pxbLayer(myid) == 1) {
1199  BRL1RX.push_back(rechit_x);
1200  BRL1RY.push_back(rechit_y);
1201  BRL1SX.push_back(sim_x);
1202  BRL1SY.push_back(sim_y);
1203  }
1204  if (tTopo->pxbLayer(myid) == 2) {
1205  BRL2RX.push_back(rechit_x);
1206  BRL2RY.push_back(rechit_y);
1207  BRL2SX.push_back(sim_x);
1208  BRL2SY.push_back(sim_y);
1209  }
1210  if (tTopo->pxbLayer(myid) == 3) {
1211  BRL3RX.push_back(rechit_x);
1212  BRL3RY.push_back(rechit_y);
1213  BRL3SX.push_back(sim_x);
1214  BRL3SY.push_back(sim_y);
1215  }
1216  }
1217 
1218  // get Forward pixels
1219  if (subid == sdPxlFwd) {
1220 
1221  ++nPxlFwd;
1222 
1223  if (tTopo->pxfDisk(myid) == 1) {
1224  if (tTopo->pxfSide(myid) == 1) {
1225  FWD1nRX.push_back(rechit_x);
1226  FWD1nRY.push_back(rechit_y);
1227  FWD1nSX.push_back(sim_x);
1228  FWD1nSY.push_back(sim_y);
1229  }
1230  if (tTopo->pxfSide(myid) == 2) {
1231  FWD1pRX.push_back(rechit_x);
1232  FWD1pRY.push_back(rechit_y);
1233  FWD1pSX.push_back(sim_x);
1234  FWD1pSY.push_back(sim_y);
1235  }
1236  }
1237  if (tTopo->pxfDisk(myid) == 2) {
1238  if (tTopo->pxfSide(myid) == 1) {
1239  FWD2nRX.push_back(rechit_x);
1240  FWD2nRY.push_back(rechit_y);
1241  FWD2nSX.push_back(sim_x);
1242  FWD2nSY.push_back(sim_y);
1243  }
1244  if (tTopo->pxfSide(myid) == 2) {
1245  FWD2pRX.push_back(rechit_x);
1246  FWD2pRY.push_back(rechit_y);
1247  FWD2pSX.push_back(sim_x);
1248  FWD2pSY.push_back(sim_y);
1249  }
1250  }
1251  }
1252  } // end matched emtpy
1253  } // <-----end rechit loop
1254  } // <------ end detunit loop
1255 
1256 
1257  if (verbosity > 1) {
1258  eventout += "\n Number of BrlPixelRecHits collected:...... ";
1259  eventout += nPxlBrl;
1260  }
1261 
1262  if (verbosity > 1) {
1263  eventout += "\n Number of FrwdPixelRecHits collected:..... ";
1264  eventout += nPxlFwd;
1265  }
1266 
1267  if (verbosity > 0)
1268  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1269 
1270  return;
1271 }
1272 
1274 {
1275  std::string MsgLoggerCat = "GlobalRecHitsProducer_storeTrk";
1276 
1277  if (verbosity > 2) {
1278 
1279  // strip output
1280  TString eventout("\n nTIBL1 = ");
1281  eventout += TIBL1RX.size();
1282  for (unsigned int i = 0; i < TIBL1RX.size(); ++i) {
1283  eventout += "\n (RX, RY, SX, SY) = (";
1284  eventout += TIBL1RX[i];
1285  eventout += ", ";
1286  eventout += TIBL1RY[i];
1287  eventout += ", ";
1288  eventout += TIBL1SX[i];
1289  eventout += ", ";
1290  eventout += TIBL1SY[i];
1291  eventout += ")";
1292  }
1293  eventout += "\n nTIBL2 = ";
1294  eventout += TIBL2RX.size();
1295  for (unsigned int i = 0; i < TIBL2RX.size(); ++i) {
1296  eventout += "\n (RX, RY, SX, SY) = (";
1297  eventout += TIBL2RX[i];
1298  eventout += ", ";
1299  eventout += TIBL2RY[i];
1300  eventout += ", ";
1301  eventout += TIBL2SX[i];
1302  eventout += ", ";
1303  eventout += TIBL2SY[i];
1304  eventout += ")";
1305  }
1306  eventout += "\n nTIBL3 = ";
1307  eventout += TIBL3RX.size();
1308  for (unsigned int i = 0; i < TIBL3RX.size(); ++i) {
1309  eventout += "\n (RX, RY, SX, SY) = (";
1310  eventout += TIBL3RX[i];
1311  eventout += ", ";
1312  eventout += TIBL3RY[i];
1313  eventout += ", ";
1314  eventout += TIBL3SX[i];
1315  eventout += ", ";
1316  eventout += TIBL3SY[i];
1317  eventout += ")";
1318  }
1319  eventout += "\n nTIBL4 = ";
1320  eventout += TIBL4RX.size();
1321  for (unsigned int i = 0; i < TIBL4RX.size(); ++i) {
1322  eventout += "\n (RX, RY, SX, SY) = (";
1323  eventout += TIBL4RX[i];
1324  eventout += ", ";
1325  eventout += TIBL4RY[i];
1326  eventout += ", ";
1327  eventout += TIBL4SX[i];
1328  eventout += ", ";
1329  eventout += TIBL4SY[i];
1330  eventout += ")";
1331  }
1332  eventout += "\n nTOBL1 = ";
1333  eventout += TOBL1RX.size();
1334  for (unsigned int i = 0; i < TOBL1RX.size(); ++i) {
1335  eventout += "\n (RX, RY, SX, SY) = (";
1336  eventout += TOBL1RX[i];
1337  eventout += ", ";
1338  eventout += TOBL1RY[i];
1339  eventout += ", ";
1340  eventout += TOBL1SX[i];
1341  eventout += ", ";
1342  eventout += TOBL1SY[i];
1343  eventout += ")";
1344  }
1345  eventout += "\n nTOBL2 = ";
1346  eventout += TOBL2RX.size();
1347  for (unsigned int i = 0; i < TOBL2RX.size(); ++i) {
1348  eventout += "\n (RX, RY, SX, SY) = (";
1349  eventout += TOBL2RX[i];
1350  eventout += ", ";
1351  eventout += TOBL2RY[i];
1352  eventout += ", ";
1353  eventout += TOBL2SX[i];
1354  eventout += ", ";
1355  eventout += TOBL2SY[i];
1356  eventout += ")";
1357  }
1358  eventout += "\n nTOBL3 = ";
1359  eventout += TOBL3RX.size();
1360  for (unsigned int i = 0; i < TOBL3RX.size(); ++i) {
1361  eventout += "\n (RX, RY, SX, SY) = (";
1362  eventout += TOBL3RX[i];
1363  eventout += ", ";
1364  eventout += TOBL3RY[i];
1365  eventout += ", ";
1366  eventout += TOBL3SX[i];
1367  eventout += ", ";
1368  eventout += TOBL3SY[i];
1369  eventout += ")";
1370  }
1371  eventout += "\n nTOBL4 = ";
1372  eventout += TOBL4RX.size();
1373  for (unsigned int i = 0; i < TOBL4RX.size(); ++i) {
1374  eventout += "\n (RX, RY, SX, SY) = (";
1375  eventout += TOBL4RX[i];
1376  eventout += ", ";
1377  eventout += TOBL4RY[i];
1378  eventout += ", ";
1379  eventout += TOBL4SX[i];
1380  eventout += ", ";
1381  eventout += TOBL4SY[i];
1382  eventout += ")";
1383  }
1384  eventout += "\n nTIDW1 = ";
1385  eventout += TIDW1RX.size();
1386  for (unsigned int i = 0; i < TIDW1RX.size(); ++i) {
1387  eventout += "\n (RX, RY, SX, SY) = (";
1388  eventout += TIDW1RX[i];
1389  eventout += ", ";
1390  eventout += TIDW1RY[i];
1391  eventout += ", ";
1392  eventout += TIDW1SX[i];
1393  eventout += ", ";
1394  eventout += TIDW1SY[i];
1395  eventout += ")";
1396  }
1397  eventout += "\n nTIDW2 = ";
1398  eventout += TIDW2RX.size();
1399  for (unsigned int i = 0; i < TIDW2RX.size(); ++i) {
1400  eventout += "\n (RX, RY, SX, SY) = (";
1401  eventout += TIDW2RX[i];
1402  eventout += ", ";
1403  eventout += TIDW2RY[i];
1404  eventout += ", ";
1405  eventout += TIDW2SX[i];
1406  eventout += ", ";
1407  eventout += TIDW2SY[i];
1408  eventout += ")";
1409  }
1410  eventout += "\n nTIDW3 = ";
1411  eventout += TIDW3RX.size();
1412  for (unsigned int i = 0; i < TIDW3RX.size(); ++i) {
1413  eventout += "\n (RX, RY, SX, SY) = (";
1414  eventout += TIDW3RX[i];
1415  eventout += ", ";
1416  eventout += TIDW3RY[i];
1417  eventout += ", ";
1418  eventout += TIDW3SX[i];
1419  eventout += ", ";
1420  eventout += TIDW3SY[i];
1421  eventout += ")";
1422  }
1423  eventout += "\n nTECW1 = ";
1424  eventout += TECW1RX.size();
1425  for (unsigned int i = 0; i < TECW1RX.size(); ++i) {
1426  eventout += "\n (RX, RY, SX, SY) = (";
1427  eventout += TECW1RX[i];
1428  eventout += ", ";
1429  eventout += TECW1RY[i];
1430  eventout += ", ";
1431  eventout += TECW1SX[i];
1432  eventout += ", ";
1433  eventout += TECW1SY[i];
1434  eventout += ")";
1435  }
1436  eventout += "\n nTECW2 = ";
1437  eventout += TECW2RX.size();
1438  for (unsigned int i = 0; i < TECW2RX.size(); ++i) {
1439  eventout += "\n (RX, RY, SX, SY) = (";
1440  eventout += TECW2RX[i];
1441  eventout += ", ";
1442  eventout += TECW2RY[i];
1443  eventout += ", ";
1444  eventout += TECW2SX[i];
1445  eventout += ", ";
1446  eventout += TECW2SY[i];
1447  eventout += ")";
1448  }
1449  eventout += "\n nTECW3 = ";
1450  eventout += TECW3RX.size();
1451  for (unsigned int i = 0; i < TECW3RX.size(); ++i) {
1452  eventout += "\n (RX, RY, SX, SY) = (";
1453  eventout += TECW3RX[i];
1454  eventout += ", ";
1455  eventout += TECW3RY[i];
1456  eventout += ", ";
1457  eventout += TECW3SX[i];
1458  eventout += ", ";
1459  eventout += TECW3SY[i];
1460  eventout += ")";
1461  }
1462  eventout += "\n nTECW4 = ";
1463  eventout += TECW4RX.size();
1464  for (unsigned int i = 0; i < TECW4RX.size(); ++i) {
1465  eventout += "\n (RX, RY, SX, SY) = (";
1466  eventout += TECW4RX[i];
1467  eventout += ", ";
1468  eventout += TECW4RY[i];
1469  eventout += ", ";
1470  eventout += TECW4SX[i];
1471  eventout += ", ";
1472  eventout += TECW4SY[i];
1473  eventout += ")";
1474  }
1475  eventout += "\n nTECW5 = ";
1476  eventout += TECW5RX.size();
1477  for (unsigned int i = 0; i < TECW5RX.size(); ++i) {
1478  eventout += "\n (RX, RY, SX, SY) = (";
1479  eventout += TECW5RX[i];
1480  eventout += ", ";
1481  eventout += TECW5RY[i];
1482  eventout += ", ";
1483  eventout += TECW5SX[i];
1484  eventout += ", ";
1485  eventout += TECW5SY[i];
1486  eventout += ")";
1487  }
1488  eventout += "\n nTECW6 = ";
1489  eventout += TECW6RX.size();
1490  for (unsigned int i = 0; i < TECW6RX.size(); ++i) {
1491  eventout += "\n (RX, RY, SX, SY) = (";
1492  eventout += TECW6RX[i];
1493  eventout += ", ";
1494  eventout += TECW6RY[i];
1495  eventout += ", ";
1496  eventout += TECW6SX[i];
1497  eventout += ", ";
1498  eventout += TECW6SY[i];
1499  eventout += ")";
1500  }
1501  eventout += "\n nTECW7 = ";
1502  eventout += TECW7RX.size();
1503  for (unsigned int i = 0; i < TECW7RX.size(); ++i) {
1504  eventout += "\n (RX, RY, SX, SY) = (";
1505  eventout += TECW7RX[i];
1506  eventout += ", ";
1507  eventout += TECW7RY[i];
1508  eventout += ", ";
1509  eventout += TECW7SX[i];
1510  eventout += ", ";
1511  eventout += TECW7SY[i];
1512  eventout += ")";
1513  }
1514  eventout += "\n nTECW8 = ";
1515  eventout += TECW8RX.size();
1516  for (unsigned int i = 0; i < TECW8RX.size(); ++i) {
1517  eventout += "\n (RX, RY, SX, SY) = (";
1518  eventout += TECW8RX[i];
1519  eventout += ", ";
1520  eventout += TECW8RY[i];
1521  eventout += ", ";
1522  eventout += TECW8SX[i];
1523  eventout += ", ";
1524  eventout += TECW8SY[i];
1525  eventout += ")";
1526  }
1527 
1528  // pixel output
1529  eventout += "\n nBRL1 = ";
1530  eventout += BRL1RX.size();
1531  for (unsigned int i = 0; i < BRL1RX.size(); ++i) {
1532  eventout += "\n (RX, RY, SX, SY) = (";
1533  eventout += BRL1RX[i];
1534  eventout += ", ";
1535  eventout += BRL1RY[i];
1536  eventout += ", ";
1537  eventout += BRL1SX[i];
1538  eventout += ", ";
1539  eventout += BRL1SY[i];
1540  eventout += ")";
1541  }
1542  eventout += "\n nBRL2 = ";
1543  eventout += BRL2RX.size();
1544  for (unsigned int i = 0; i < BRL2RX.size(); ++i) {
1545  eventout += "\n (RX, RY, SX, SY) = (";
1546  eventout += BRL2RX[i];
1547  eventout += ", ";
1548  eventout += BRL2RY[i];
1549  eventout += ", ";
1550  eventout += BRL2SX[i];
1551  eventout += ", ";
1552  eventout += BRL2SY[i];
1553  eventout += ")";
1554  }
1555  eventout += "\n nBRL3 = ";
1556  eventout += BRL3RX.size();
1557  for (unsigned int i = 0; i < BRL3RX.size(); ++i) {
1558  eventout += "\n (RX, RY, SX, SY) = (";
1559  eventout += BRL3RX[i];
1560  eventout += ", ";
1561  eventout += BRL3RY[i];
1562  eventout += ", ";
1563  eventout += BRL3SX[i];
1564  eventout += ", ";
1565  eventout += BRL3SY[i];
1566  eventout += ")";
1567  }
1568  eventout += "\n nFWD1p = ";
1569  eventout += FWD1pRX.size();
1570  for (unsigned int i = 0; i < FWD1pRX.size(); ++i) {
1571  eventout += "\n (RX, RY, SX, SY) = (";
1572  eventout += FWD1pRX[i];
1573  eventout += ", ";
1574  eventout += FWD1pRY[i];
1575  eventout += ", ";
1576  eventout += FWD1pSX[i];
1577  eventout += ", ";
1578  eventout += FWD1pSY[i];
1579  eventout += ")";
1580  }
1581  eventout += "\n nFWD1n = ";
1582  eventout += FWD1nRX.size();
1583  for (unsigned int i = 0; i < FWD1nRX.size(); ++i) {
1584  eventout += "\n (RX, RY, SX, SY) = (";
1585  eventout += FWD1nRX[i];
1586  eventout += ", ";
1587  eventout += FWD1nRY[i];
1588  eventout += ", ";
1589  eventout += FWD1nSX[i];
1590  eventout += ", ";
1591  eventout += FWD1nSY[i];
1592  eventout += ")";
1593  }
1594  eventout += "\n nFWD2p = ";
1595  eventout += FWD2pRX.size();
1596  for (unsigned int i = 0; i < FWD2pRX.size(); ++i) {
1597  eventout += "\n (RX, RY, SX, SY) = (";
1598  eventout += FWD2pRX[i];
1599  eventout += ", ";
1600  eventout += FWD2pRY[i];
1601  eventout += ", ";
1602  eventout += FWD2pSX[i];
1603  eventout += ", ";
1604  eventout += FWD2pSY[i];
1605  eventout += ")";
1606  }
1607  eventout += "\n nFWD2p = ";
1608  eventout += FWD2nRX.size();
1609  for (unsigned int i = 0; i < FWD2nRX.size(); ++i) {
1610  eventout += "\n (RX, RY, SX, SY) = (";
1611  eventout += FWD2nRX[i];
1612  eventout += ", ";
1613  eventout += FWD2nRY[i];
1614  eventout += ", ";
1615  eventout += FWD2nSX[i];
1616  eventout += ", ";
1617  eventout += FWD2nSY[i];
1618  eventout += ")";
1619  }
1620 
1621  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1622  }
1623 
1624  // strip output
1644 
1645  // pixel output
1653 
1654  return;
1655 }
1656 
1658  const edm::EventSetup& iSetup)
1659 {
1660  std::string MsgLoggerCat = "GlobalRecHitsProducer_fillMuon";
1661 
1662  TString eventout;
1663  if (verbosity > 0)
1664  eventout = "\nGathering info:";
1665 
1666  // get DT information
1668  iSetup.get<MuonGeometryRecord>().get(dtGeom);
1669  if (!dtGeom.isValid()) {
1670  edm::LogWarning(MsgLoggerCat)
1671  << "Unable to find DTMuonGeometryRecord in event!";
1672  return;
1673  }
1674 
1676  iEvent.getByLabel(MuDTSimSrc_, dtsimHits);
1677  if (!dtsimHits.isValid()) {
1678  edm::LogWarning(MsgLoggerCat)
1679  << "Unable to find dtsimHits in event!";
1680  return;
1681  }
1682 
1683  std::map<DTWireId, edm::PSimHitContainer> simHitsPerWire =
1685 
1687  iEvent.getByLabel(MuDTSrc_, dtRecHits);
1688  if (!dtRecHits.isValid()) {
1689  edm::LogWarning(MsgLoggerCat)
1690  << "Unable to find dtRecHits in event!";
1691  return;
1692  }
1693 
1694  std::map<DTWireId, std::vector<DTRecHit1DPair> > recHitsPerWire =
1695  map1DRecHitsPerWire(dtRecHits.product());
1696 
1697 
1698  int nDt = compute(dtGeom.product(), simHitsPerWire, recHitsPerWire, 1);
1699 
1700  if (verbosity > 1) {
1701  eventout += "\n Number of DtMuonRecHits collected:........ ";
1702  eventout += nDt;
1703  }
1704 
1705  // get CSC Strip information
1706  // get map of sim hits
1707  theMap.clear();
1708  //edm::Handle<CrossingFrame> cf;
1710  //iEvent.getByType(cf);
1711  //if (!cf.isValid()) {
1712  // edm::LogWarning(MsgLoggerCat)
1713  // << "Unable to find CrossingFrame in event!";
1714  // return;
1715  //}
1716  //MixCollection<PSimHit> simHits(cf.product(), "MuonCSCHits");
1717  iEvent.getByLabel("mix","MuonCSCHits",cf);
1718  if (!cf.isValid()) {
1719  edm::LogWarning(MsgLoggerCat)
1720  << "Unable to find muo CSC crossingFrame in event!";
1721  return;
1722  }
1724 
1725  // arrange the hits by detUnit
1726  for(MixCollection<PSimHit>::MixItr hitItr = simHits.begin();
1727  hitItr != simHits.end(); ++hitItr)
1728  {
1729  theMap[hitItr->detUnitId()].push_back(*hitItr);
1730  }
1731 
1732  // get geometry
1734  iSetup.get<MuonGeometryRecord>().get(hGeom);
1735  if (!hGeom.isValid()) {
1736  edm::LogWarning(MsgLoggerCat)
1737  << "Unable to find CSCMuonGeometryRecord in event!";
1738  return;
1739  }
1740  const CSCGeometry *theCSCGeometry = &*hGeom;
1741 
1742  // get rechits
1744  iEvent.getByLabel(MuCSCSrc_, hRecHits);
1745  if (!hRecHits.isValid()) {
1746  edm::LogWarning(MsgLoggerCat)
1747  << "Unable to find CSC RecHits in event!";
1748  return;
1749  }
1750  const CSCRecHit2DCollection *cscRecHits = hRecHits.product();
1751 
1752  int nCSC = 0;
1753  for (CSCRecHit2DCollection::const_iterator recHitItr = cscRecHits->begin();
1754  recHitItr != cscRecHits->end(); ++recHitItr) {
1755 
1756  int detId = (*recHitItr).cscDetId().rawId();
1757 
1759  std::map<int, edm::PSimHitContainer>::const_iterator mapItr =
1760  theMap.find(detId);
1761  if (mapItr != theMap.end()) {
1762  simHits = mapItr->second;
1763  }
1764 
1765  if (simHits.size() == 1) {
1766  ++nCSC;
1767 
1768  const GeomDetUnit* detUnit =
1769  theCSCGeometry->idToDetUnit(CSCDetId(detId));
1770  const CSCLayer *layer = dynamic_cast<const CSCLayer *>(detUnit);
1771 
1772  int chamberType = layer->chamber()->specs()->chamberType();
1773  plotResolution(simHits[0], *recHitItr, layer, chamberType);
1774  }
1775  }
1776 
1777  if (verbosity > 1) {
1778  eventout += "\n Number of CSCRecHits collected:........... ";
1779  eventout += nCSC;
1780  }
1781 
1782  // get RPC information
1783  std::map<double, int> mapsim, maprec;
1784  std::map<int, double> nmapsim, nmaprec;
1785 
1787  iSetup.get<MuonGeometryRecord>().get(rpcGeom);
1788  if (!rpcGeom.isValid()) {
1789  edm::LogWarning(MsgLoggerCat)
1790  << "Unable to find RPCMuonGeometryRecord in event!";
1791  return;
1792  }
1793 
1795  iEvent.getByLabel(MuRPCSimSrc_, simHit);
1796  if (!simHit.isValid()) {
1797  edm::LogWarning(MsgLoggerCat)
1798  << "Unable to find RPCSimHit in event!";
1799  return;
1800  }
1801 
1803  iEvent.getByLabel(MuRPCSrc_, recHit);
1804  if (!simHit.isValid()) {
1805  edm::LogWarning(MsgLoggerCat)
1806  << "Unable to find RPCRecHit in event!";
1807  return;
1808  }
1809 
1810  int nRPC = 0;
1812  int nrec = 0;
1813  for (recIt = recHit->begin(); recIt != recHit->end(); ++recIt) {
1814  RPCDetId Rid = (RPCDetId)(*recIt).rpcId();
1815  const RPCRoll *roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(Rid));
1816  if (roll->isForward()) {
1817 
1818  if (verbosity > 1) {
1819  eventout += "\n Number of RPCRecHits collected:........... ";
1820  eventout += nRPC;
1821  }
1822 
1823  if (verbosity > 0)
1824  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1825  return;
1826  }
1827  nrec = nrec + 1;
1828  LocalPoint rhitlocal = (*recIt).localPosition();
1829  double rhitlocalx = rhitlocal.x();
1830  maprec[rhitlocalx] = nrec;
1831  }
1832 
1833  int i = 0;
1834  for (std::map<double,int>::iterator iter = maprec.begin();
1835  iter != maprec.end(); ++iter) {
1836  i = i + 1;
1837  nmaprec[i] = (*iter).first;
1838  }
1839 
1840  edm::PSimHitContainer::const_iterator simIt;
1841  int nsim = 0;
1842  for (simIt = simHit->begin(); simIt != simHit->end(); simIt++) {
1843  int ptype = (*simIt).particleType();
1844  //RPCDetId Rsid = (RPCDetId)(*simIt).detUnitId();
1845  if (ptype == 13 || ptype == -13) {
1846  nsim = nsim + 1;
1847  LocalPoint shitlocal = (*simIt).localPosition();
1848  double shitlocalx = shitlocal.x();
1849  mapsim[shitlocalx] = nsim;
1850  }
1851  }
1852 
1853  i = 0;
1854  for (std::map<double,int>::iterator iter = mapsim.begin();
1855  iter != mapsim.end(); ++iter) {
1856  i = i + 1;
1857  nmapsim[i] = (*iter).first;
1858  }
1859 
1860  if (nsim == nrec) {
1861  for (int r = 0; r < nsim; r++) {
1862  ++nRPC;
1863  RPCRHX.push_back(nmaprec[r+1]);
1864  RPCSHX.push_back(nmapsim[r+1]);
1865  }
1866  }
1867 
1868  if (verbosity > 1) {
1869  eventout += "\n Number of RPCRecHits collected:........... ";
1870  eventout += nRPC;
1871  }
1872 
1873  if (verbosity > 0)
1874  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1875 
1876  return;
1877 }
1878 
1880 {
1881  std::string MsgLoggerCat = "GlobalRecHitsProducer_storeMuon";
1882 
1883  if (verbosity > 2) {
1884 
1885  // dt output
1886  TString eventout("\n nDT = ");
1887  eventout += DTRHD.size();
1888  for (unsigned int i = 0; i < DTRHD.size(); ++i) {
1889  eventout += "\n (RHD, SHD) = (";
1890  eventout += DTRHD[i];
1891  eventout += ", ";
1892  eventout += DTSHD[i];
1893  eventout += ")";
1894  }
1895 
1896  // CSC Strip
1897  eventout += "\n nCSC = ";
1898  eventout += CSCRHPHI.size();
1899  for (unsigned int i = 0; i < CSCRHPHI.size(); ++i) {
1900  eventout += "\n (rhphi, rhperp, shphi) = (";
1901  eventout += CSCRHPHI[i];
1902  eventout += ", ";
1903  eventout += CSCRHPERP[i];
1904  eventout += ", ";
1905  eventout += CSCSHPHI[i];
1906  eventout += ")";
1907  }
1908 
1909  // RPC
1910  eventout += "\n nRPC = ";
1911  eventout += RPCRHX.size();
1912  for (unsigned int i = 0; i < RPCRHX.size(); ++i) {
1913  eventout += "\n (rhx, shx) = (";
1914  eventout += RPCRHX[i];
1915  eventout += ", ";
1916  eventout += RPCSHX[i];
1917  eventout += ")";
1918  }
1919 
1920  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1921  }
1922 
1923  product.putDTRecHits(DTRHD,DTSHD);
1924 
1926 
1927  product.putRPCRecHits(RPCRHX,RPCSHX);
1928 
1929  return;
1930 }
1931 
1933 {
1934  std::string MsgLoggerCat = "GlobalRecHitsProducer_clear";
1935 
1936  if (verbosity > 0)
1937  edm::LogInfo(MsgLoggerCat)
1938  << "Clearing event holders";
1939 
1940  // reset electromagnetic info
1941  // EE info
1942  EERE.clear();
1943  EESHE.clear();
1944  // EB info
1945  EBRE.clear();
1946  EBSHE.clear();
1947  // ES info
1948  ESRE.clear();
1949  ESSHE.clear();
1950 
1951  // reset HCal Info
1952  HBCalREC.clear();
1953  HBCalR.clear();
1954  HBCalSHE.clear();
1955  HECalREC.clear();
1956  HECalR.clear();
1957  HECalSHE.clear();
1958  HOCalREC.clear();
1959  HOCalR.clear();
1960  HOCalSHE.clear();
1961  HFCalREC.clear();
1962  HFCalR.clear();
1963  HFCalSHE.clear();
1964 
1965  // reset Track Info
1966  TIBL1RX.clear();
1967  TIBL2RX.clear();
1968  TIBL3RX.clear();
1969  TIBL4RX.clear();
1970  TIBL1RY.clear();
1971  TIBL2RY.clear();
1972  TIBL3RY.clear();
1973  TIBL4RY.clear();
1974  TIBL1SX.clear();
1975  TIBL2SX.clear();
1976  TIBL3SX.clear();
1977  TIBL4SX.clear();
1978  TIBL1SY.clear();
1979  TIBL2SY.clear();
1980  TIBL3SY.clear();
1981  TIBL4SY.clear();
1982 
1983  TOBL1RX.clear();
1984  TOBL2RX.clear();
1985  TOBL3RX.clear();
1986  TOBL4RX.clear();
1987  TOBL1RY.clear();
1988  TOBL2RY.clear();
1989  TOBL3RY.clear();
1990  TOBL4RY.clear();
1991  TOBL1SX.clear();
1992  TOBL2SX.clear();
1993  TOBL3SX.clear();
1994  TOBL4SX.clear();
1995  TOBL1SY.clear();
1996  TOBL2SY.clear();
1997  TOBL3SY.clear();
1998  TOBL4SY.clear();
1999 
2000  TIDW1RX.clear();
2001  TIDW2RX.clear();
2002  TIDW3RX.clear();
2003  TIDW1RY.clear();
2004  TIDW2RY.clear();
2005  TIDW3RY.clear();
2006  TIDW1SX.clear();
2007  TIDW2SX.clear();
2008  TIDW3SX.clear();
2009  TIDW1SY.clear();
2010  TIDW2SY.clear();
2011  TIDW3SY.clear();
2012 
2013  TECW1RX.clear();
2014  TECW2RX.clear();
2015  TECW3RX.clear();
2016  TECW4RX.clear();
2017  TECW5RX.clear();
2018  TECW6RX.clear();
2019  TECW7RX.clear();
2020  TECW8RX.clear();
2021  TECW1RY.clear();
2022  TECW2RY.clear();
2023  TECW3RY.clear();
2024  TECW4RY.clear();
2025  TECW5RY.clear();
2026  TECW6RY.clear();
2027  TECW7RY.clear();
2028  TECW8RY.clear();
2029  TECW1SX.clear();
2030  TECW2SX.clear();
2031  TECW3SX.clear();
2032  TECW4SX.clear();
2033  TECW5SX.clear();
2034  TECW6SX.clear();
2035  TECW7SX.clear();
2036  TECW8SX.clear();
2037  TECW1SY.clear();
2038  TECW2SY.clear();
2039  TECW3SY.clear();
2040  TECW4SY.clear();
2041  TECW5SY.clear();
2042  TECW6SY.clear();
2043  TECW7SY.clear();
2044  TECW8SY.clear();
2045 
2046  BRL1RX.clear();
2047  BRL1RY.clear();
2048  BRL1SX.clear();
2049  BRL1SY.clear();
2050  BRL2RX.clear();
2051  BRL2RY.clear();
2052  BRL2SX.clear();
2053  BRL2SY.clear();
2054  BRL3RX.clear();
2055  BRL3RY.clear();
2056  BRL3SX.clear();
2057  BRL3SY.clear();
2058 
2059  FWD1pRX.clear();
2060  FWD1pRY.clear();
2061  FWD1pSX.clear();
2062  FWD1pSY.clear();
2063  FWD1nRX.clear();
2064  FWD1nRY.clear();
2065  FWD1nSX.clear();
2066  FWD1nSY.clear();
2067  FWD2pRX.clear();
2068  FWD2pRY.clear();
2069  FWD2pSX.clear();
2070  FWD2pSY.clear();
2071  FWD2nRX.clear();
2072  FWD2nRY.clear();
2073  FWD2nSX.clear();
2074  FWD2nSY.clear();
2075 
2076  //muon clear
2077  DTRHD.clear();
2078  DTSHD.clear();
2079 
2080  CSCRHPHI.clear();
2081  CSCRHPERP.clear();
2082  CSCSHPHI.clear();
2083 
2084  RPCRHX.clear();
2085  RPCSHX.clear();
2086 
2087  return;
2088 }
2089 
2090 //needed by to do the residual for matched hits in SiStrip
2091 std::pair<LocalPoint,LocalVector>
2093  const StripGeomDetUnit* stripDet,
2094  const BoundPlane& plane)
2095 {
2096 
2097  const StripTopology& topol = stripDet->specificTopology();
2098  GlobalPoint globalpos= stripDet->surface().toGlobal(hit.localPosition());
2099  LocalPoint localHit = plane.toLocal(globalpos);
2100  //track direction
2101  LocalVector locdir=hit.localDirection();
2102  //rotate track in new frame
2103 
2104  GlobalVector globaldir= stripDet->surface().toGlobal(locdir);
2105  LocalVector dir=plane.toLocal(globaldir);
2106  float scale = -localHit.z() / dir.z();
2107 
2108  LocalPoint projectedPos = localHit + scale*dir;
2109 
2110  float selfAngle = topol.stripAngle( topol.strip( hit.localPosition()));
2111 
2112  // vector along strip in hit frame
2113  LocalVector stripDir( sin(selfAngle), cos(selfAngle), 0);
2114 
2115  LocalVector
2116  localStripDir(plane.toLocal(stripDet->surface().toGlobal(stripDir)));
2117 
2118  return std::pair<LocalPoint,LocalVector>( projectedPos, localStripDir);
2119 }
2120 
2121 // Return a map between DTRecHit1DPair and wireId
2122 std::map<DTWireId, std::vector<DTRecHit1DPair> >
2124  dt1DRecHitPairs) {
2125  std::map<DTWireId, std::vector<DTRecHit1DPair> > ret;
2126 
2127  for(DTRecHitCollection::const_iterator rechit = dt1DRecHitPairs->begin();
2128  rechit != dt1DRecHitPairs->end(); rechit++) {
2129  ret[(*rechit).wireId()].push_back(*rechit);
2130  }
2131 
2132  return ret;
2133 }
2134 
2135 // Compute SimHit distance from wire (cm)
2137  DTWireId wireId,
2138  const PSimHit& hit) {
2139  float xwire = layer->specificTopology().wirePosition(wireId.wire());
2140  LocalPoint entryP = hit.entryPoint();
2141  LocalPoint exitP = hit.exitPoint();
2142  float xEntry = entryP.x()-xwire;
2143  float xExit = exitP.x()-xwire;
2144 
2145  //FIXME: check...
2146  return fabs(xEntry - (entryP.z()*(xExit-xEntry))/(exitP.z()-entryP.z()));
2147 }
2148 
2149 // Find the RecHit closest to the muon SimHit
2150 template <typename type>
2151 const type*
2153  DTWireId wireId,
2154  const std::vector<type>& recHits,
2155  const float simHitDist) {
2156  float res = 99999;
2157  const type* theBestRecHit = 0;
2158  // Loop over RecHits within the cell
2159  for(typename std::vector<type>::const_iterator recHit = recHits.begin();
2160  recHit != recHits.end();
2161  recHit++) {
2162  float distTmp = recHitDistFromWire(*recHit, layer);
2163  if(fabs(distTmp-simHitDist) < res) {
2164  res = fabs(distTmp-simHitDist);
2165  theBestRecHit = &(*recHit);
2166  }
2167  } // End of loop over RecHits within the cell
2168 
2169  return theBestRecHit;
2170 }
2171 
2172 // Compute the distance from wire (cm) of a hits in a DTRecHit1DPair
2173 float
2175  const DTLayer* layer) {
2176  // Compute the rechit distance from wire
2177  return fabs(hitPair.localPosition(DTEnums::Left).x() -
2178  hitPair.localPosition(DTEnums::Right).x())/2.;
2179 }
2180 
2181 // Compute the distance from wire (cm) of a hits in a DTRecHit1D
2182 float
2184  const DTLayer* layer) {
2185  return fabs(recHit.localPosition().x() -
2186  layer->specificTopology().wirePosition(recHit.wireId().wire()));
2187 }
2188 
2189 template <typename type>
2191  const std::map<DTWireId, std::vector<PSimHit> >&
2192  _simHitsPerWire,
2193  const std::map<DTWireId, std::vector<type> >&
2194  _recHitsPerWire,
2195  int step) {
2196 
2197  std::map<DTWireId, std::vector<PSimHit> > simHitsPerWire = _simHitsPerWire;
2198  std::map<DTWireId, std::vector<type> > recHitsPerWire = _recHitsPerWire;
2199  int nDt = 0;
2200  // Loop over cells with a muon SimHit
2201  for(std::map<DTWireId, std::vector<PSimHit> >::const_iterator wireAndSHits =
2202  simHitsPerWire.begin();
2203  wireAndSHits != simHitsPerWire.end();
2204  wireAndSHits++) {
2205  DTWireId wireId = (*wireAndSHits).first;
2206  std::vector<PSimHit> simHitsInCell = (*wireAndSHits).second;
2207 
2208  // Get the layer
2209  const DTLayer* layer = dtGeom->layer(wireId);
2210 
2211  // Look for a mu hit in the cell
2212  const PSimHit* muSimHit = DTHitQualityUtils::findMuSimHit(simHitsInCell);
2213  if (muSimHit==0) {
2214  continue; // Skip this cell
2215  }
2216 
2217  // Find the distance of the simhit from the wire
2218  float simHitWireDist = simHitDistFromWire(layer, wireId, *muSimHit);
2219  // Skip simhits out of the cell
2220  if(simHitWireDist>2.1) {
2221  continue; // Skip this cell
2222  }
2223  //GlobalPoint simHitGlobalPos = layer->toGlobal(muSimHit->localPosition());
2224 
2225  // Look for RecHits in the same cell
2226  if(recHitsPerWire.find(wireId) == recHitsPerWire.end()) {
2227  continue; // No RecHit found in this cell
2228  } else {
2229 
2230  // vector<type> recHits = (*wireAndRecHits).second;
2231  std::vector<type> recHits = recHitsPerWire[wireId];
2232 
2233  // Find the best RecHit
2234  const type* theBestRecHit =
2235  findBestRecHit(layer, wireId, recHits, simHitWireDist);
2236 
2237  float recHitWireDist = recHitDistFromWire(*theBestRecHit, layer);
2238 
2239  ++nDt;
2240 
2241  DTRHD.push_back(recHitWireDist);
2242  DTSHD.push_back(simHitWireDist);
2243 
2244  } // find rechits
2245  } // loop over simhits
2246 
2247  return nDt;
2248 }
2249 
2250 void
2252  const CSCRecHit2D & recHit,
2253  const CSCLayer * layer,
2254  int chamberType) {
2255  GlobalPoint simHitPos = layer->toGlobal(simHit.localPosition());
2256  GlobalPoint recHitPos = layer->toGlobal(recHit.localPosition());
2257 
2258  CSCRHPHI.push_back(recHitPos.phi());
2259  CSCRHPERP.push_back(recHitPos.perp());
2260  CSCSHPHI.push_back(simHitPos.phi());
2261 }
2262 
2263 //define this as a plug-in
2264 //DEFINE_FWK_MODULE(GlobalRecHitsProducer);
void putTIBL1RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
RunNumber_t run() const
Definition: EventID.h:42
GlobalRecHitsProducer(const edm::ParameterSet &)
void getManyByType(std::vector< Handle< PROD > > &results) const
Definition: Event.h:395
type
Definition: HCALResponse.h:21
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:114
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:44
T getUntrackedParameter(std::string const &, T const &) const
float wirePosition(int wireNumber) const
Returns the x position in the layer of a given wire number.
Definition: DTTopology.cc:88
void getAllProvenance(std::vector< Provenance const * > &provenances) const
Definition: Event.cc:78
int i
Definition: DBlmapReader.cc:9
const type * findBestRecHit(const DTLayer *layer, DTWireId wireId, const std::vector< type > &recHits, const float simHitDist)
float simHitDistFromWire(const DTLayer *layer, DTWireId wireId, const PSimHit &hit)
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
void putHFCalRecHits(const std::vector< float > &rec, const std::vector< float > &r, const std::vector< float > &she)
std::vector< PCaloHit > PCaloHitContainer
T perp() const
Definition: PV3DBase.h:72
int ihf
unsigned int tibLayer(const DetId &id) const
void putHECalRecHits(const std::vector< float > &rec, const std::vector< float > &r, const std::vector< float > &she)
virtual float stripAngle(float strip) const =0
void putFWD1nRecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
void putTOBL4RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
static const int sdHcalOut
#define PI
float recHitDistFromWire(const DTRecHit1DPair &hitPair, const DTLayer *layer)
void putHBCalRecHits(const std::vector< float > &rec, const std::vector< float > &r, const std::vector< float > &she)
std::pair< LocalPoint, LocalVector > projectHit(const PSimHit &hit, const StripGeomDetUnit *stripDet, const BoundPlane &plane)
std::vector< PSimHit > matched
void putTECW6RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
unsigned int pxfDisk(const DetId &id) const
void putBRL1RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
void putTECW8RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
void fillHCal(edm::Event &, const edm::EventSetup &)
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:47
void putEECalRecHits(const std::vector< float > &re, const std::vector< float > &she)
static const int sdSiTID
void putTECW5RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
void putTECW7RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
void putCSCRecHits(const std::vector< float > &rhphi, const std::vector< float > &rhperp, const std::vector< float > &shphi)
void putTECW2RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
std::vector< T >::const_iterator const_iterator
T y() const
Definition: PV3DBase.h:63
static const PSimHit * findMuSimHit(const edm::PSimHitContainer &hits)
Select the SimHit from a muon in a vector of SimHits.
void putTIDW1RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
void putRPCRecHits(const std::vector< float > &rhx, const std::vector< float > &shx)
void putTECW1RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
unsigned int tidWheel(const DetId &id) const
void storeMuon(PGlobalRecHit &)
virtual const StripTopology & specificTopology() const
Returns a reference to the strip proxy topology.
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:35
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
void storeECal(PGlobalRecHit &)
std::map< int, edm::PSimHitContainer > theMap
dictionary map
Definition: Association.py:205
void putBRL2RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
virtual float strip(const LocalPoint &) const =0
virtual void produce(edm::Event &, const edm::EventSetup &)
const DTLayer * layer(DTLayerId id) const
Return a layer given its id.
Definition: DTGeometry.cc:112
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
C::const_iterator const_iterator
constant access iterator type
Definition: RangeMap.h:45
const DTTopology & specificTopology() const
Definition: DTLayer.cc:44
int iEvent
Definition: GenABIO.cc:243
Local3DPoint exitPoint() const
Exit point in the local Det frame.
Definition: PSimHit.h:38
void putTOBL1RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
void putTIDW2RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
void putFWD2nRecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
static const int sdSiTIB
const CSCChamberSpecs * specs() const
Definition: CSCChamber.h:42
static const int sdPxlBrl
Local3DPoint localPosition() const
Definition: PSimHit.h:44
virtual const GeomDetUnit * idToDetUnit(DetId) const
Return the pointer to the GeomDetUnit corresponding to a given DetId.
Definition: CSCGeometry.cc:93
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:94
void putHOCalRecHits(const std::vector< float > &rec, const std::vector< float > &r, const std::vector< float > &she)
void putTECW3RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
void putTOBL3RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
T sqrt(T t)
Definition: SSEVec.h:48
void putTIBL2RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
T z() const
Definition: PV3DBase.h:64
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
std::map< uint32_t, float, std::less< uint32_t > > MapType
void plotResolution(const PSimHit &simHit, const CSCRecHit2D &recHit, const CSCLayer *layer, int chamberType)
static const int sdSiTOB
void putEBCalRecHits(const std::vector< float > &re, const std::vector< float > &she)
void putBRL3RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
bool isValid() const
Definition: HandleBase.h:76
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:39
virtual const GeomDet * idToDet(DetId) const
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:361
void storeTrk(PGlobalRecHit &)
virtual LocalPoint localPosition() const
Return the 3-dimensional local position.
Definition: DTRecHit1D.h:62
unsigned int pxbLayer(const DetId &id) const
const_iterator end() const
void putTIBL4RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
LocalVector localDirection() const
Obsolete. Same as momentumAtEntry().unit(), for backward compatibility.
Definition: PSimHit.h:52
int wire() const
Return the wire number.
Definition: DTWireId.h:58
void putESCalRecHits(const std::vector< float > &re, const std::vector< float > &she)
Definition: DetId.h:20
void putFWD1pRecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
void fillMuon(edm::Event &, const edm::EventSetup &)
int chamberType() const
void putDTRecHits(const std::vector< float > &rhd, const std::vector< float > &shd)
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:62
static const int sdHcalFwd
tuple simHits
Definition: trackerHits.py:16
T const * product() const
Definition: Handle.h:74
static const int sdHcalBrl
std::string const & label() const
Definition: InputTag.h:42
T eta() const
Definition: PV3DBase.h:76
void storeHCal(PGlobalRecHit &)
void putTOBL2RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
ESHandle< TrackerGeometry > geometry
edm::EventID id() const
Definition: EventBase.h:56
static const int sdSiTEC
iterator find(key_type k)
iterator end()
Definition: DetSetNew.h:59
std::vector< PSimHit > associateHit(const TrackingRecHit &thit)
unsigned int pxfSide(const DetId &id) const
std::vector< PSimHit > PSimHitContainer
static const int sdPxlFwd
dbl *** dir
Definition: mlp_gen.cc:35
static std::map< DTWireId, edm::PSimHitContainer > mapSimHitsPerWire(const edm::PSimHitContainer &simhits)
bool isValid() const
Definition: ESHandle.h:37
void fillTrk(edm::Event &, const edm::EventSetup &)
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
bool isForward() const
Definition: RPCRoll.cc:98
T x() const
Definition: PV3DBase.h:62
unsigned int tecWheel(const DetId &id) const
int compute(const DTGeometry *dtGeom, const std::map< DTWireId, std::vector< PSimHit > > &simHitsPerWire, const std::map< DTWireId, std::vector< type > > &recHitsPerWire, int step)
std::string const & instance() const
Definition: InputTag.h:43
Local3DPoint entryPoint() const
Entry point in the local Det frame.
Definition: PSimHit.h:35
void putTECW4RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
static const int sdHcalEC
const CSCChamber * chamber() const
Definition: CSCLayer.h:52
void fillECal(edm::Event &, const edm::EventSetup &)
void putFWD2pRecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
const GeomDetUnit * stereoDet() const
Definition: GluedGeomDet.h:21
const_iterator begin() const
void putTIDW3RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
unsigned int tobLayer(const DetId &id) const
void putTIBL3RecHits(const std::vector< float > &rx, const std::vector< float > &ry, const std::vector< float > &sx, const std::vector< float > &sy)
std::map< DTWireId, std::vector< DTRecHit1DPair > > map1DRecHitsPerWire(const DTRecHitCollection *dt1DRecHitPairs)
virtual LocalPoint localPosition() const
DTWireId wireId() const
Return the wireId.
Definition: DTRecHit1D.h:109
iterator begin()
Definition: DetSetNew.h:56