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