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 maxHOEnergy = 0.;
550  Double_t maxHFEnergy = 0.;
551 
552  Double_t maxHBPhi = -1000.;
553  Double_t maxHEPhi = -1000.;
554  Double_t maxHOPhi = -1000.;
555  Double_t maxHFPhi = -1000.;
556 
557  Double_t maxHBEta = -1000.;
558  Double_t maxHEEta = -1000.;
559  Double_t maxHOEta = -1000.;
560  Double_t maxHFEta = -1000.;
561 
562  Double_t PI = 3.141592653589;
563 
565  // get HBHE information
567  std::vector<edm::Handle<HBHERecHitCollection> > hbhe;
568  iEvent.getManyByType(hbhe);
569  if (!hbhe[0].isValid()) {
570  edm::LogWarning(MsgLoggerCat)
571  << "Unable to find any HBHERecHitCollections in event!";
572  return;
573  }
574  std::vector<edm::Handle<HBHERecHitCollection> >::iterator ihbhe;
575 
576  int iHB = 0;
577  int iHE = 0;
578  for (ihbhe = hbhe.begin(); ihbhe != hbhe.end(); ++ihbhe) {
579 
580  // find max values
581  for (HBHERecHitCollection::const_iterator jhbhe = (*ihbhe)->begin();
582  jhbhe != (*ihbhe)->end(); ++jhbhe) {
583 
584  HcalDetId cell(jhbhe->id());
585 
586  if (cell.subdet() == sdHcalBrl) {
587 
588  const CaloCellGeometry* cellGeometry =
589  geometry->getSubdetectorGeometry (cell)->getGeometry (cell) ;
590  double fEta = cellGeometry->getPosition().eta () ;
591  double fPhi = cellGeometry->getPosition().phi () ;
592  if ( (jhbhe->energy()) > maxHBEnergy ) {
593  maxHBEnergy = jhbhe->energy();
594  maxHOEnergy = maxHBEnergy;
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  std::string MsgLoggerCat = "GlobalRecHitsProducer_fillTrk";
859 
860  TString eventout;
861  if (verbosity > 0)
862  eventout = "\nGathering info:";
863 
864  // get strip information
866  iEvent.getByLabel(SiStripSrc_, rechitsmatched);
867  if (!rechitsmatched.isValid()) {
868  edm::LogWarning(MsgLoggerCat)
869  << "Unable to find stripmatchedrechits in event!";
870  return;
871  }
872 
873  TrackerHitAssociator associate(iEvent,conf_);
874 
876  iSetup.get<TrackerDigiGeometryRecord>().get(pDD);
877  if (!pDD.isValid()) {
878  edm::LogWarning(MsgLoggerCat)
879  << "Unable to find TrackerDigiGeometry in event!";
880  return;
881  }
882  const TrackerGeometry &tracker(*pDD);
883 
884  int nStripBrl = 0, nStripFwd = 0;
885 
886  // loop over det units
887  for (TrackerGeometry::DetContainer::const_iterator it = pDD->dets().begin();
888  it != pDD->dets().end(); ++it) {
889 
890  uint32_t myid = ((*it)->geographicalId()).rawId();
891  DetId detid = ((*it)->geographicalId());
892 
893  //loop over rechits-matched in the same subdetector
894  SiStripMatchedRecHit2DCollection::const_iterator rechitmatchedMatch = rechitsmatched->find(detid);
895 
896  if (rechitmatchedMatch != rechitsmatched->end()) {
897  SiStripMatchedRecHit2DCollection::DetSet rechitmatchedRange = *rechitmatchedMatch;
898  SiStripMatchedRecHit2DCollection::DetSet::const_iterator rechitmatchedRangeIteratorBegin = rechitmatchedRange.begin();
899  SiStripMatchedRecHit2DCollection::DetSet::const_iterator rechitmatchedRangeIteratorEnd = rechitmatchedRange.end();
900  SiStripMatchedRecHit2DCollection::DetSet::const_iterator itermatched = rechitmatchedRangeIteratorBegin;
901 
902  for ( itermatched = rechitmatchedRangeIteratorBegin;
903  itermatched != rechitmatchedRangeIteratorEnd;
904  ++itermatched) {
905 
906  SiStripMatchedRecHit2D const rechit = *itermatched;
907  LocalPoint position = rechit.localPosition();
908 
909  float mindist = 999999.;
910  float distx = 999999.;
911  float disty = 999999.;
912  float dist = 999999.;
913  std::pair<LocalPoint,LocalVector> closestPair;
914  matched.clear();
915 
916  float rechitmatchedx = position.x();
917  float rechitmatchedy = position.y();
918 
919  matched = associate.associateHit(rechit);
920 
921  if (!matched.empty()) {
922  //project simhit;
923  const GluedGeomDet* gluedDet =
924  (const GluedGeomDet*)tracker.idToDet(rechit.geographicalId());
925  const StripGeomDetUnit* partnerstripdet =
926  (StripGeomDetUnit*) gluedDet->stereoDet();
927  std::pair<LocalPoint,LocalVector> hitPair;
928 
929  for(std::vector<PSimHit>::const_iterator m = matched.begin();
930  m != matched.end(); m++){
931  //project simhit;
932  hitPair = projectHit((*m),partnerstripdet,gluedDet->surface());
933  distx = fabs(rechitmatchedx - hitPair.first.x());
934  disty = fabs(rechitmatchedy - hitPair.first.y());
935  dist = sqrt(distx*distx+disty*disty);
936 
937  if(dist < mindist){
938  mindist = dist;
939  closestPair = hitPair;
940  }
941  }
942 
943  // get TIB
944  if (detid.subdetId() == sdSiTIB) {
945 
946  TIBDetId tibid(myid);
947  ++nStripBrl;
948 
949  if (tibid.layer() == 1) {
950  TIBL1RX.push_back(rechitmatchedx);
951  TIBL1RY.push_back(rechitmatchedy);
952  TIBL1SX.push_back(closestPair.first.x());
953  TIBL1SY.push_back(closestPair.first.y());
954  }
955  if (tibid.layer() == 2) {
956  TIBL2RX.push_back(rechitmatchedx);
957  TIBL2RY.push_back(rechitmatchedy);
958  TIBL2SX.push_back(closestPair.first.x());
959  TIBL2SY.push_back(closestPair.first.y());
960  }
961  if (tibid.layer() == 3) {
962  TIBL3RX.push_back(rechitmatchedx);
963  TIBL3RY.push_back(rechitmatchedy);
964  TIBL3SX.push_back(closestPair.first.x());
965  TIBL3SY.push_back(closestPair.first.y());
966  }
967  if (tibid.layer() == 4) {
968  TIBL4RX.push_back(rechitmatchedx);
969  TIBL4RY.push_back(rechitmatchedy);
970  TIBL4SX.push_back(closestPair.first.x());
971  TIBL4SY.push_back(closestPair.first.y());
972  }
973  }
974 
975  // get TOB
976  if (detid.subdetId() == sdSiTOB) {
977 
978  TOBDetId tobid(myid);
979  ++nStripBrl;
980 
981  if (tobid.layer() == 1) {
982  TOBL1RX.push_back(rechitmatchedx);
983  TOBL1RY.push_back(rechitmatchedy);
984  TOBL1SX.push_back(closestPair.first.x());
985  TOBL1SY.push_back(closestPair.first.y());
986  }
987  if (tobid.layer() == 2) {
988  TOBL2RX.push_back(rechitmatchedx);
989  TOBL2RY.push_back(rechitmatchedy);
990  TOBL2SX.push_back(closestPair.first.x());
991  TOBL2SY.push_back(closestPair.first.y());
992  }
993  if (tobid.layer() == 3) {
994  TOBL3RX.push_back(rechitmatchedx);
995  TOBL3RY.push_back(rechitmatchedy);
996  TOBL3SX.push_back(closestPair.first.x());
997  TOBL3SY.push_back(closestPair.first.y());
998  }
999  if (tobid.layer() == 4) {
1000  TOBL4RX.push_back(rechitmatchedx);
1001  TOBL4RY.push_back(rechitmatchedy);
1002  TOBL4SX.push_back(closestPair.first.x());
1003  TOBL4SY.push_back(closestPair.first.y());
1004  }
1005  }
1006 
1007  // get TID
1008  if (detid.subdetId() == sdSiTID) {
1009 
1010  TIDDetId tidid(myid);
1011  ++nStripFwd;
1012 
1013  if (tidid.wheel() == 1) {
1014  TIDW1RX.push_back(rechitmatchedx);
1015  TIDW1RY.push_back(rechitmatchedy);
1016  TIDW1SX.push_back(closestPair.first.x());
1017  TIDW1SY.push_back(closestPair.first.y());
1018  }
1019  if (tidid.wheel() == 2) {
1020  TIDW2RX.push_back(rechitmatchedx);
1021  TIDW2RY.push_back(rechitmatchedy);
1022  TIDW2SX.push_back(closestPair.first.x());
1023  TIDW2SY.push_back(closestPair.first.y());
1024  }
1025  if (tidid.wheel() == 3) {
1026  TIDW3RX.push_back(rechitmatchedx);
1027  TIDW3RY.push_back(rechitmatchedy);
1028  TIDW3SX.push_back(closestPair.first.x());
1029  TIDW3SY.push_back(closestPair.first.y());
1030  }
1031  }
1032 
1033  // get TEC
1034  if (detid.subdetId() == sdSiTEC) {
1035 
1036  TECDetId tecid(myid);
1037  ++nStripFwd;
1038 
1039  if (tecid.wheel() == 1) {
1040  TECW1RX.push_back(rechitmatchedx);
1041  TECW1RY.push_back(rechitmatchedy);
1042  TECW1SX.push_back(closestPair.first.x());
1043  TECW1SY.push_back(closestPair.first.y());
1044  }
1045  if (tecid.wheel() == 2) {
1046  TECW2RX.push_back(rechitmatchedx);
1047  TECW2RY.push_back(rechitmatchedy);
1048  TECW2SX.push_back(closestPair.first.x());
1049  TECW2SY.push_back(closestPair.first.y());
1050  }
1051  if (tecid.wheel() == 3) {
1052  TECW3RX.push_back(rechitmatchedx);
1053  TECW3RY.push_back(rechitmatchedy);
1054  TECW3SX.push_back(closestPair.first.x());
1055  TECW3SY.push_back(closestPair.first.y());
1056  }
1057  if (tecid.wheel() == 4) {
1058  TECW4RX.push_back(rechitmatchedx);
1059  TECW4RY.push_back(rechitmatchedy);
1060  TECW4SX.push_back(closestPair.first.x());
1061  TECW4SY.push_back(closestPair.first.y());
1062  }
1063  if (tecid.wheel() == 5) {
1064  TECW5RX.push_back(rechitmatchedx);
1065  TECW5RY.push_back(rechitmatchedy);
1066  TECW5SX.push_back(closestPair.first.x());
1067  TECW5SY.push_back(closestPair.first.y());
1068  }
1069  if (tecid.wheel() == 6) {
1070  TECW6RX.push_back(rechitmatchedx);
1071  TECW6RY.push_back(rechitmatchedy);
1072  TECW6SX.push_back(closestPair.first.x());
1073  TECW6SY.push_back(closestPair.first.y());
1074  }
1075  if (tecid.wheel() == 7) {
1076  TECW7RX.push_back(rechitmatchedx);
1077  TECW7RY.push_back(rechitmatchedy);
1078  TECW7SX.push_back(closestPair.first.x());
1079  TECW7SY.push_back(closestPair.first.y());
1080  }
1081  if (tecid.wheel() == 8) {
1082  TECW8RX.push_back(rechitmatchedx);
1083  TECW8RY.push_back(rechitmatchedy);
1084  TECW8SX.push_back(closestPair.first.x());
1085  TECW8SY.push_back(closestPair.first.y());
1086  }
1087  }
1088 
1089  } // end if matched empty
1090  }
1091  }
1092  } // end loop over det units
1093 
1094  if (verbosity > 1) {
1095  eventout += "\n Number of BrlStripRecHits collected:...... ";
1096  eventout += nStripBrl;
1097  }
1098 
1099  if (verbosity > 1) {
1100  eventout += "\n Number of FrwdStripRecHits collected:..... ";
1101  eventout += nStripFwd;
1102  }
1103 
1104  // get pixel information
1105  //Get RecHits
1107  iEvent.getByLabel(SiPxlSrc_, recHitColl);
1108  if (!recHitColl.isValid()) {
1109  edm::LogWarning(MsgLoggerCat)
1110  << "Unable to find SiPixelRecHitCollection in event!";
1111  return;
1112  }
1113 
1114  //Get event setup
1116  iSetup.get<TrackerDigiGeometryRecord>().get(geom);
1117  if (!geom.isValid()) {
1118  edm::LogWarning(MsgLoggerCat)
1119  << "Unable to find TrackerDigiGeometry in event!";
1120  return;
1121  }
1122  //const TrackerGeometry& theTracker(*geom);
1123 
1124  int nPxlBrl = 0, nPxlFwd = 0;
1125  //iterate over detunits
1126  for (TrackerGeometry::DetContainer::const_iterator it = geom->dets().begin();
1127  it != geom->dets().end(); ++it) {
1128 
1129  uint32_t myid = ((*it)->geographicalId()).rawId();
1130  DetId detId = ((*it)->geographicalId());
1131  int subid = detId.subdetId();
1132 
1133  if (! ((subid == sdPxlBrl) || (subid == sdPxlFwd))) continue;
1134 
1135  //const PixelGeomDetUnit * theGeomDet =
1136  // dynamic_cast<const PixelGeomDetUnit*>(theTracker.idToDet(detId) );
1137 
1138  SiPixelRecHitCollection::const_iterator pixeldet = recHitColl->find(detId);
1139  if (pixeldet == recHitColl->end()) continue;
1140  SiPixelRecHitCollection::DetSet pixelrechitRange = *pixeldet;
1141  SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorBegin = pixelrechitRange.begin();
1142  SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorEnd = pixelrechitRange.end();
1143  SiPixelRecHitCollection::DetSet::const_iterator pixeliter = pixelrechitRangeIteratorBegin;
1144  std::vector<PSimHit> matched;
1145 
1146  //----Loop over rechits for this detId
1147  for ( ; pixeliter != pixelrechitRangeIteratorEnd; ++pixeliter) {
1148 
1149  matched.clear();
1150  matched = associate.associateHit(*pixeliter);
1151 
1152  if ( !matched.empty() ) {
1153 
1154  float closest = 9999.9;
1155  //std::vector<PSimHit>::const_iterator closestit = matched.begin();
1156  LocalPoint lp = pixeliter->localPosition();
1157  float rechit_x = lp.x();
1158  float rechit_y = lp.y();
1159 
1160  float sim_x = 0.;
1161  float sim_y = 0.;
1162 
1163  //loop over sim hits and fill closet
1164  for (std::vector<PSimHit>::const_iterator m = matched.begin();
1165  m != matched.end(); ++m) {
1166 
1167  float sim_x1 = (*m).entryPoint().x();
1168  float sim_x2 = (*m).exitPoint().x();
1169  float sim_xpos = 0.5*(sim_x1+sim_x2);
1170 
1171  float sim_y1 = (*m).entryPoint().y();
1172  float sim_y2 = (*m).exitPoint().y();
1173  float sim_ypos = 0.5*(sim_y1+sim_y2);
1174 
1175  float x_res = fabs(sim_xpos - rechit_x);
1176  float y_res = fabs(sim_ypos - rechit_y);
1177 
1178  float dist = sqrt(x_res*x_res + y_res*y_res);
1179 
1180  if ( dist < closest ) {
1181  closest = dist;
1182  sim_x = sim_xpos;
1183  sim_y = sim_ypos;
1184  }
1185  } // end sim hit loop
1186 
1187  // get Barrel pixels
1188  if (subid == sdPxlBrl) {
1189  PXBDetId bdetid(myid);
1190  ++nPxlBrl;
1191 
1192  if (bdetid.layer() == 1) {
1193  BRL1RX.push_back(rechit_x);
1194  BRL1RY.push_back(rechit_y);
1195  BRL1SX.push_back(sim_x);
1196  BRL1SY.push_back(sim_y);
1197  }
1198  if (bdetid.layer() == 2) {
1199  BRL2RX.push_back(rechit_x);
1200  BRL2RY.push_back(rechit_y);
1201  BRL2SX.push_back(sim_x);
1202  BRL2SY.push_back(sim_y);
1203  }
1204  if (bdetid.layer() == 3) {
1205  BRL3RX.push_back(rechit_x);
1206  BRL3RY.push_back(rechit_y);
1207  BRL3SX.push_back(sim_x);
1208  BRL3SY.push_back(sim_y);
1209  }
1210  }
1211 
1212  // get Forward pixels
1213  if (subid == sdPxlFwd) {
1214  PXFDetId fdetid(myid);
1215  ++nPxlFwd;
1216 
1217  if (fdetid.disk() == 1) {
1218  if (fdetid.side() == 1) {
1219  FWD1nRX.push_back(rechit_x);
1220  FWD1nRY.push_back(rechit_y);
1221  FWD1nSX.push_back(sim_x);
1222  FWD1nSY.push_back(sim_y);
1223  }
1224  if (fdetid.side() == 2) {
1225  FWD1pRX.push_back(rechit_x);
1226  FWD1pRY.push_back(rechit_y);
1227  FWD1pSX.push_back(sim_x);
1228  FWD1pSY.push_back(sim_y);
1229  }
1230  }
1231  if (fdetid.disk() == 2) {
1232  if (fdetid.side() == 1) {
1233  FWD2nRX.push_back(rechit_x);
1234  FWD2nRY.push_back(rechit_y);
1235  FWD2nSX.push_back(sim_x);
1236  FWD2nSY.push_back(sim_y);
1237  }
1238  if (fdetid.side() == 2) {
1239  FWD2pRX.push_back(rechit_x);
1240  FWD2pRY.push_back(rechit_y);
1241  FWD2pSX.push_back(sim_x);
1242  FWD2pSY.push_back(sim_y);
1243  }
1244  }
1245  }
1246  } // end matched emtpy
1247  } // <-----end rechit loop
1248  } // <------ end detunit loop
1249 
1250 
1251  if (verbosity > 1) {
1252  eventout += "\n Number of BrlPixelRecHits collected:...... ";
1253  eventout += nPxlBrl;
1254  }
1255 
1256  if (verbosity > 1) {
1257  eventout += "\n Number of FrwdPixelRecHits collected:..... ";
1258  eventout += nPxlFwd;
1259  }
1260 
1261  if (verbosity > 0)
1262  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1263 
1264  return;
1265 }
1266 
1268 {
1269  std::string MsgLoggerCat = "GlobalRecHitsProducer_storeTrk";
1270 
1271  if (verbosity > 2) {
1272 
1273  // strip output
1274  TString eventout("\n nTIBL1 = ");
1275  eventout += TIBL1RX.size();
1276  for (unsigned int i = 0; i < TIBL1RX.size(); ++i) {
1277  eventout += "\n (RX, RY, SX, SY) = (";
1278  eventout += TIBL1RX[i];
1279  eventout += ", ";
1280  eventout += TIBL1RY[i];
1281  eventout += ", ";
1282  eventout += TIBL1SX[i];
1283  eventout += ", ";
1284  eventout += TIBL1SY[i];
1285  eventout += ")";
1286  }
1287  eventout += "\n nTIBL2 = ";
1288  eventout += TIBL2RX.size();
1289  for (unsigned int i = 0; i < TIBL2RX.size(); ++i) {
1290  eventout += "\n (RX, RY, SX, SY) = (";
1291  eventout += TIBL2RX[i];
1292  eventout += ", ";
1293  eventout += TIBL2RY[i];
1294  eventout += ", ";
1295  eventout += TIBL2SX[i];
1296  eventout += ", ";
1297  eventout += TIBL2SY[i];
1298  eventout += ")";
1299  }
1300  eventout += "\n nTIBL3 = ";
1301  eventout += TIBL3RX.size();
1302  for (unsigned int i = 0; i < TIBL3RX.size(); ++i) {
1303  eventout += "\n (RX, RY, SX, SY) = (";
1304  eventout += TIBL3RX[i];
1305  eventout += ", ";
1306  eventout += TIBL3RY[i];
1307  eventout += ", ";
1308  eventout += TIBL3SX[i];
1309  eventout += ", ";
1310  eventout += TIBL3SY[i];
1311  eventout += ")";
1312  }
1313  eventout += "\n nTIBL4 = ";
1314  eventout += TIBL4RX.size();
1315  for (unsigned int i = 0; i < TIBL4RX.size(); ++i) {
1316  eventout += "\n (RX, RY, SX, SY) = (";
1317  eventout += TIBL4RX[i];
1318  eventout += ", ";
1319  eventout += TIBL4RY[i];
1320  eventout += ", ";
1321  eventout += TIBL4SX[i];
1322  eventout += ", ";
1323  eventout += TIBL4SY[i];
1324  eventout += ")";
1325  }
1326  eventout += "\n nTOBL1 = ";
1327  eventout += TOBL1RX.size();
1328  for (unsigned int i = 0; i < TOBL1RX.size(); ++i) {
1329  eventout += "\n (RX, RY, SX, SY) = (";
1330  eventout += TOBL1RX[i];
1331  eventout += ", ";
1332  eventout += TOBL1RY[i];
1333  eventout += ", ";
1334  eventout += TOBL1SX[i];
1335  eventout += ", ";
1336  eventout += TOBL1SY[i];
1337  eventout += ")";
1338  }
1339  eventout += "\n nTOBL2 = ";
1340  eventout += TOBL2RX.size();
1341  for (unsigned int i = 0; i < TOBL2RX.size(); ++i) {
1342  eventout += "\n (RX, RY, SX, SY) = (";
1343  eventout += TOBL2RX[i];
1344  eventout += ", ";
1345  eventout += TOBL2RY[i];
1346  eventout += ", ";
1347  eventout += TOBL2SX[i];
1348  eventout += ", ";
1349  eventout += TOBL2SY[i];
1350  eventout += ")";
1351  }
1352  eventout += "\n nTOBL3 = ";
1353  eventout += TOBL3RX.size();
1354  for (unsigned int i = 0; i < TOBL3RX.size(); ++i) {
1355  eventout += "\n (RX, RY, SX, SY) = (";
1356  eventout += TOBL3RX[i];
1357  eventout += ", ";
1358  eventout += TOBL3RY[i];
1359  eventout += ", ";
1360  eventout += TOBL3SX[i];
1361  eventout += ", ";
1362  eventout += TOBL3SY[i];
1363  eventout += ")";
1364  }
1365  eventout += "\n nTOBL4 = ";
1366  eventout += TOBL4RX.size();
1367  for (unsigned int i = 0; i < TOBL4RX.size(); ++i) {
1368  eventout += "\n (RX, RY, SX, SY) = (";
1369  eventout += TOBL4RX[i];
1370  eventout += ", ";
1371  eventout += TOBL4RY[i];
1372  eventout += ", ";
1373  eventout += TOBL4SX[i];
1374  eventout += ", ";
1375  eventout += TOBL4SY[i];
1376  eventout += ")";
1377  }
1378  eventout += "\n nTIDW1 = ";
1379  eventout += TIDW1RX.size();
1380  for (unsigned int i = 0; i < TIDW1RX.size(); ++i) {
1381  eventout += "\n (RX, RY, SX, SY) = (";
1382  eventout += TIDW1RX[i];
1383  eventout += ", ";
1384  eventout += TIDW1RY[i];
1385  eventout += ", ";
1386  eventout += TIDW1SX[i];
1387  eventout += ", ";
1388  eventout += TIDW1SY[i];
1389  eventout += ")";
1390  }
1391  eventout += "\n nTIDW2 = ";
1392  eventout += TIDW2RX.size();
1393  for (unsigned int i = 0; i < TIDW2RX.size(); ++i) {
1394  eventout += "\n (RX, RY, SX, SY) = (";
1395  eventout += TIDW2RX[i];
1396  eventout += ", ";
1397  eventout += TIDW2RY[i];
1398  eventout += ", ";
1399  eventout += TIDW2SX[i];
1400  eventout += ", ";
1401  eventout += TIDW2SY[i];
1402  eventout += ")";
1403  }
1404  eventout += "\n nTIDW3 = ";
1405  eventout += TIDW3RX.size();
1406  for (unsigned int i = 0; i < TIDW3RX.size(); ++i) {
1407  eventout += "\n (RX, RY, SX, SY) = (";
1408  eventout += TIDW3RX[i];
1409  eventout += ", ";
1410  eventout += TIDW3RY[i];
1411  eventout += ", ";
1412  eventout += TIDW3SX[i];
1413  eventout += ", ";
1414  eventout += TIDW3SY[i];
1415  eventout += ")";
1416  }
1417  eventout += "\n nTECW1 = ";
1418  eventout += TECW1RX.size();
1419  for (unsigned int i = 0; i < TECW1RX.size(); ++i) {
1420  eventout += "\n (RX, RY, SX, SY) = (";
1421  eventout += TECW1RX[i];
1422  eventout += ", ";
1423  eventout += TECW1RY[i];
1424  eventout += ", ";
1425  eventout += TECW1SX[i];
1426  eventout += ", ";
1427  eventout += TECW1SY[i];
1428  eventout += ")";
1429  }
1430  eventout += "\n nTECW2 = ";
1431  eventout += TECW2RX.size();
1432  for (unsigned int i = 0; i < TECW2RX.size(); ++i) {
1433  eventout += "\n (RX, RY, SX, SY) = (";
1434  eventout += TECW2RX[i];
1435  eventout += ", ";
1436  eventout += TECW2RY[i];
1437  eventout += ", ";
1438  eventout += TECW2SX[i];
1439  eventout += ", ";
1440  eventout += TECW2SY[i];
1441  eventout += ")";
1442  }
1443  eventout += "\n nTECW3 = ";
1444  eventout += TECW3RX.size();
1445  for (unsigned int i = 0; i < TECW3RX.size(); ++i) {
1446  eventout += "\n (RX, RY, SX, SY) = (";
1447  eventout += TECW3RX[i];
1448  eventout += ", ";
1449  eventout += TECW3RY[i];
1450  eventout += ", ";
1451  eventout += TECW3SX[i];
1452  eventout += ", ";
1453  eventout += TECW3SY[i];
1454  eventout += ")";
1455  }
1456  eventout += "\n nTECW4 = ";
1457  eventout += TECW4RX.size();
1458  for (unsigned int i = 0; i < TECW4RX.size(); ++i) {
1459  eventout += "\n (RX, RY, SX, SY) = (";
1460  eventout += TECW4RX[i];
1461  eventout += ", ";
1462  eventout += TECW4RY[i];
1463  eventout += ", ";
1464  eventout += TECW4SX[i];
1465  eventout += ", ";
1466  eventout += TECW4SY[i];
1467  eventout += ")";
1468  }
1469  eventout += "\n nTECW5 = ";
1470  eventout += TECW5RX.size();
1471  for (unsigned int i = 0; i < TECW5RX.size(); ++i) {
1472  eventout += "\n (RX, RY, SX, SY) = (";
1473  eventout += TECW5RX[i];
1474  eventout += ", ";
1475  eventout += TECW5RY[i];
1476  eventout += ", ";
1477  eventout += TECW5SX[i];
1478  eventout += ", ";
1479  eventout += TECW5SY[i];
1480  eventout += ")";
1481  }
1482  eventout += "\n nTECW6 = ";
1483  eventout += TECW6RX.size();
1484  for (unsigned int i = 0; i < TECW6RX.size(); ++i) {
1485  eventout += "\n (RX, RY, SX, SY) = (";
1486  eventout += TECW6RX[i];
1487  eventout += ", ";
1488  eventout += TECW6RY[i];
1489  eventout += ", ";
1490  eventout += TECW6SX[i];
1491  eventout += ", ";
1492  eventout += TECW6SY[i];
1493  eventout += ")";
1494  }
1495  eventout += "\n nTECW7 = ";
1496  eventout += TECW7RX.size();
1497  for (unsigned int i = 0; i < TECW7RX.size(); ++i) {
1498  eventout += "\n (RX, RY, SX, SY) = (";
1499  eventout += TECW7RX[i];
1500  eventout += ", ";
1501  eventout += TECW7RY[i];
1502  eventout += ", ";
1503  eventout += TECW7SX[i];
1504  eventout += ", ";
1505  eventout += TECW7SY[i];
1506  eventout += ")";
1507  }
1508  eventout += "\n nTECW8 = ";
1509  eventout += TECW8RX.size();
1510  for (unsigned int i = 0; i < TECW8RX.size(); ++i) {
1511  eventout += "\n (RX, RY, SX, SY) = (";
1512  eventout += TECW8RX[i];
1513  eventout += ", ";
1514  eventout += TECW8RY[i];
1515  eventout += ", ";
1516  eventout += TECW8SX[i];
1517  eventout += ", ";
1518  eventout += TECW8SY[i];
1519  eventout += ")";
1520  }
1521 
1522  // pixel output
1523  eventout += "\n nBRL1 = ";
1524  eventout += BRL1RX.size();
1525  for (unsigned int i = 0; i < BRL1RX.size(); ++i) {
1526  eventout += "\n (RX, RY, SX, SY) = (";
1527  eventout += BRL1RX[i];
1528  eventout += ", ";
1529  eventout += BRL1RY[i];
1530  eventout += ", ";
1531  eventout += BRL1SX[i];
1532  eventout += ", ";
1533  eventout += BRL1SY[i];
1534  eventout += ")";
1535  }
1536  eventout += "\n nBRL2 = ";
1537  eventout += BRL2RX.size();
1538  for (unsigned int i = 0; i < BRL2RX.size(); ++i) {
1539  eventout += "\n (RX, RY, SX, SY) = (";
1540  eventout += BRL2RX[i];
1541  eventout += ", ";
1542  eventout += BRL2RY[i];
1543  eventout += ", ";
1544  eventout += BRL2SX[i];
1545  eventout += ", ";
1546  eventout += BRL2SY[i];
1547  eventout += ")";
1548  }
1549  eventout += "\n nBRL3 = ";
1550  eventout += BRL3RX.size();
1551  for (unsigned int i = 0; i < BRL3RX.size(); ++i) {
1552  eventout += "\n (RX, RY, SX, SY) = (";
1553  eventout += BRL3RX[i];
1554  eventout += ", ";
1555  eventout += BRL3RY[i];
1556  eventout += ", ";
1557  eventout += BRL3SX[i];
1558  eventout += ", ";
1559  eventout += BRL3SY[i];
1560  eventout += ")";
1561  }
1562  eventout += "\n nFWD1p = ";
1563  eventout += FWD1pRX.size();
1564  for (unsigned int i = 0; i < FWD1pRX.size(); ++i) {
1565  eventout += "\n (RX, RY, SX, SY) = (";
1566  eventout += FWD1pRX[i];
1567  eventout += ", ";
1568  eventout += FWD1pRY[i];
1569  eventout += ", ";
1570  eventout += FWD1pSX[i];
1571  eventout += ", ";
1572  eventout += FWD1pSY[i];
1573  eventout += ")";
1574  }
1575  eventout += "\n nFWD1n = ";
1576  eventout += FWD1nRX.size();
1577  for (unsigned int i = 0; i < FWD1nRX.size(); ++i) {
1578  eventout += "\n (RX, RY, SX, SY) = (";
1579  eventout += FWD1nRX[i];
1580  eventout += ", ";
1581  eventout += FWD1nRY[i];
1582  eventout += ", ";
1583  eventout += FWD1nSX[i];
1584  eventout += ", ";
1585  eventout += FWD1nSY[i];
1586  eventout += ")";
1587  }
1588  eventout += "\n nFWD2p = ";
1589  eventout += FWD2pRX.size();
1590  for (unsigned int i = 0; i < FWD2pRX.size(); ++i) {
1591  eventout += "\n (RX, RY, SX, SY) = (";
1592  eventout += FWD2pRX[i];
1593  eventout += ", ";
1594  eventout += FWD2pRY[i];
1595  eventout += ", ";
1596  eventout += FWD2pSX[i];
1597  eventout += ", ";
1598  eventout += FWD2pSY[i];
1599  eventout += ")";
1600  }
1601  eventout += "\n nFWD2p = ";
1602  eventout += FWD2nRX.size();
1603  for (unsigned int i = 0; i < FWD2nRX.size(); ++i) {
1604  eventout += "\n (RX, RY, SX, SY) = (";
1605  eventout += FWD2nRX[i];
1606  eventout += ", ";
1607  eventout += FWD2nRY[i];
1608  eventout += ", ";
1609  eventout += FWD2nSX[i];
1610  eventout += ", ";
1611  eventout += FWD2nSY[i];
1612  eventout += ")";
1613  }
1614 
1615  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1616  }
1617 
1618  // strip output
1638 
1639  // pixel output
1647 
1648  return;
1649 }
1650 
1652  const edm::EventSetup& iSetup)
1653 {
1654  std::string MsgLoggerCat = "GlobalRecHitsProducer_fillMuon";
1655 
1656  TString eventout;
1657  if (verbosity > 0)
1658  eventout = "\nGathering info:";
1659 
1660  // get DT information
1662  iSetup.get<MuonGeometryRecord>().get(dtGeom);
1663  if (!dtGeom.isValid()) {
1664  edm::LogWarning(MsgLoggerCat)
1665  << "Unable to find DTMuonGeometryRecord in event!";
1666  return;
1667  }
1668 
1670  iEvent.getByLabel(MuDTSimSrc_, dtsimHits);
1671  if (!dtsimHits.isValid()) {
1672  edm::LogWarning(MsgLoggerCat)
1673  << "Unable to find dtsimHits in event!";
1674  return;
1675  }
1676 
1677  std::map<DTWireId, edm::PSimHitContainer> simHitsPerWire =
1679 
1681  iEvent.getByLabel(MuDTSrc_, dtRecHits);
1682  if (!dtRecHits.isValid()) {
1683  edm::LogWarning(MsgLoggerCat)
1684  << "Unable to find dtRecHits in event!";
1685  return;
1686  }
1687 
1688  std::map<DTWireId, std::vector<DTRecHit1DPair> > recHitsPerWire =
1689  map1DRecHitsPerWire(dtRecHits.product());
1690 
1691 
1692  int nDt = compute(dtGeom.product(), simHitsPerWire, recHitsPerWire, 1);
1693 
1694  if (verbosity > 1) {
1695  eventout += "\n Number of DtMuonRecHits collected:........ ";
1696  eventout += nDt;
1697  }
1698 
1699  // get CSC Strip information
1700  // get map of sim hits
1701  theMap.clear();
1702  //edm::Handle<CrossingFrame> cf;
1704  //iEvent.getByType(cf);
1705  //if (!cf.isValid()) {
1706  // edm::LogWarning(MsgLoggerCat)
1707  // << "Unable to find CrossingFrame in event!";
1708  // return;
1709  //}
1710  //MixCollection<PSimHit> simHits(cf.product(), "MuonCSCHits");
1711  iEvent.getByLabel("mix","MuonCSCHits",cf);
1712  if (!cf.isValid()) {
1713  edm::LogWarning(MsgLoggerCat)
1714  << "Unable to find muo CSC crossingFrame in event!";
1715  return;
1716  }
1718 
1719  // arrange the hits by detUnit
1720  for(MixCollection<PSimHit>::MixItr hitItr = simHits.begin();
1721  hitItr != simHits.end(); ++hitItr)
1722  {
1723  theMap[hitItr->detUnitId()].push_back(*hitItr);
1724  }
1725 
1726  // get geometry
1728  iSetup.get<MuonGeometryRecord>().get(hGeom);
1729  if (!hGeom.isValid()) {
1730  edm::LogWarning(MsgLoggerCat)
1731  << "Unable to find CSCMuonGeometryRecord in event!";
1732  return;
1733  }
1734  const CSCGeometry *theCSCGeometry = &*hGeom;
1735 
1736  // get rechits
1738  iEvent.getByLabel(MuCSCSrc_, hRecHits);
1739  if (!hRecHits.isValid()) {
1740  edm::LogWarning(MsgLoggerCat)
1741  << "Unable to find CSC RecHits in event!";
1742  return;
1743  }
1744  const CSCRecHit2DCollection *cscRecHits = hRecHits.product();
1745 
1746  int nCSC = 0;
1747  for (CSCRecHit2DCollection::const_iterator recHitItr = cscRecHits->begin();
1748  recHitItr != cscRecHits->end(); ++recHitItr) {
1749 
1750  int detId = (*recHitItr).cscDetId().rawId();
1751 
1753  std::map<int, edm::PSimHitContainer>::const_iterator mapItr =
1754  theMap.find(detId);
1755  if (mapItr != theMap.end()) {
1756  simHits = mapItr->second;
1757  }
1758 
1759  if (simHits.size() == 1) {
1760  ++nCSC;
1761 
1762  const GeomDetUnit* detUnit =
1763  theCSCGeometry->idToDetUnit(CSCDetId(detId));
1764  const CSCLayer *layer = dynamic_cast<const CSCLayer *>(detUnit);
1765 
1766  int chamberType = layer->chamber()->specs()->chamberType();
1767  plotResolution(simHits[0], *recHitItr, layer, chamberType);
1768  }
1769  }
1770 
1771  if (verbosity > 1) {
1772  eventout += "\n Number of CSCRecHits collected:........... ";
1773  eventout += nCSC;
1774  }
1775 
1776  // get RPC information
1777  std::map<double, int> mapsim, maprec;
1778  std::map<int, double> nmapsim, nmaprec;
1779 
1781  iSetup.get<MuonGeometryRecord>().get(rpcGeom);
1782  if (!rpcGeom.isValid()) {
1783  edm::LogWarning(MsgLoggerCat)
1784  << "Unable to find RPCMuonGeometryRecord in event!";
1785  return;
1786  }
1787 
1789  iEvent.getByLabel(MuRPCSimSrc_, simHit);
1790  if (!simHit.isValid()) {
1791  edm::LogWarning(MsgLoggerCat)
1792  << "Unable to find RPCSimHit in event!";
1793  return;
1794  }
1795 
1797  iEvent.getByLabel(MuRPCSrc_, recHit);
1798  if (!simHit.isValid()) {
1799  edm::LogWarning(MsgLoggerCat)
1800  << "Unable to find RPCRecHit in event!";
1801  return;
1802  }
1803 
1804  int nRPC = 0;
1806  int nrec = 0;
1807  for (recIt = recHit->begin(); recIt != recHit->end(); ++recIt) {
1808  RPCDetId Rid = (RPCDetId)(*recIt).rpcId();
1809  const RPCRoll *roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(Rid));
1810  if (roll->isForward()) {
1811 
1812  if (verbosity > 1) {
1813  eventout += "\n Number of RPCRecHits collected:........... ";
1814  eventout += nRPC;
1815  }
1816 
1817  if (verbosity > 0)
1818  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1819  return;
1820  }
1821  nrec = nrec + 1;
1822  LocalPoint rhitlocal = (*recIt).localPosition();
1823  double rhitlocalx = rhitlocal.x();
1824  maprec[rhitlocalx] = nrec;
1825  }
1826 
1827  int i = 0;
1828  for (std::map<double,int>::iterator iter = maprec.begin();
1829  iter != maprec.end(); ++iter) {
1830  i = i + 1;
1831  nmaprec[i] = (*iter).first;
1832  }
1833 
1834  edm::PSimHitContainer::const_iterator simIt;
1835  int nsim = 0;
1836  for (simIt = simHit->begin(); simIt != simHit->end(); simIt++) {
1837  int ptype = (*simIt).particleType();
1838  //RPCDetId Rsid = (RPCDetId)(*simIt).detUnitId();
1839  if (ptype == 13 || ptype == -13) {
1840  nsim = nsim + 1;
1841  LocalPoint shitlocal = (*simIt).localPosition();
1842  double shitlocalx = shitlocal.x();
1843  mapsim[shitlocalx] = nsim;
1844  }
1845  }
1846 
1847  i = 0;
1848  for (std::map<double,int>::iterator iter = mapsim.begin();
1849  iter != mapsim.end(); ++iter) {
1850  i = i + 1;
1851  nmapsim[i] = (*iter).first;
1852  }
1853 
1854  if (nsim == nrec) {
1855  for (int r = 0; r < nsim; r++) {
1856  ++nRPC;
1857  RPCRHX.push_back(nmaprec[r+1]);
1858  RPCSHX.push_back(nmapsim[r+1]);
1859  }
1860  }
1861 
1862  if (verbosity > 1) {
1863  eventout += "\n Number of RPCRecHits collected:........... ";
1864  eventout += nRPC;
1865  }
1866 
1867  if (verbosity > 0)
1868  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1869 
1870  return;
1871 }
1872 
1874 {
1875  std::string MsgLoggerCat = "GlobalRecHitsProducer_storeMuon";
1876 
1877  if (verbosity > 2) {
1878 
1879  // dt output
1880  TString eventout("\n nDT = ");
1881  eventout += DTRHD.size();
1882  for (unsigned int i = 0; i < DTRHD.size(); ++i) {
1883  eventout += "\n (RHD, SHD) = (";
1884  eventout += DTRHD[i];
1885  eventout += ", ";
1886  eventout += DTSHD[i];
1887  eventout += ")";
1888  }
1889 
1890  // CSC Strip
1891  eventout += "\n nCSC = ";
1892  eventout += CSCRHPHI.size();
1893  for (unsigned int i = 0; i < CSCRHPHI.size(); ++i) {
1894  eventout += "\n (rhphi, rhperp, shphi) = (";
1895  eventout += CSCRHPHI[i];
1896  eventout += ", ";
1897  eventout += CSCRHPERP[i];
1898  eventout += ", ";
1899  eventout += CSCSHPHI[i];
1900  eventout += ")";
1901  }
1902 
1903  // RPC
1904  eventout += "\n nRPC = ";
1905  eventout += RPCRHX.size();
1906  for (unsigned int i = 0; i < RPCRHX.size(); ++i) {
1907  eventout += "\n (rhx, shx) = (";
1908  eventout += RPCRHX[i];
1909  eventout += ", ";
1910  eventout += RPCSHX[i];
1911  eventout += ")";
1912  }
1913 
1914  edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1915  }
1916 
1917  product.putDTRecHits(DTRHD,DTSHD);
1918 
1920 
1921  product.putRPCRecHits(RPCRHX,RPCSHX);
1922 
1923  return;
1924 }
1925 
1927 {
1928  std::string MsgLoggerCat = "GlobalRecHitsProducer_clear";
1929 
1930  if (verbosity > 0)
1931  edm::LogInfo(MsgLoggerCat)
1932  << "Clearing event holders";
1933 
1934  // reset electromagnetic info
1935  // EE info
1936  EERE.clear();
1937  EESHE.clear();
1938  // EB info
1939  EBRE.clear();
1940  EBSHE.clear();
1941  // ES info
1942  ESRE.clear();
1943  ESSHE.clear();
1944 
1945  // reset HCal Info
1946  HBCalREC.clear();
1947  HBCalR.clear();
1948  HBCalSHE.clear();
1949  HECalREC.clear();
1950  HECalR.clear();
1951  HECalSHE.clear();
1952  HOCalREC.clear();
1953  HOCalR.clear();
1954  HOCalSHE.clear();
1955  HFCalREC.clear();
1956  HFCalR.clear();
1957  HFCalSHE.clear();
1958 
1959  // reset Track Info
1960  TIBL1RX.clear();
1961  TIBL2RX.clear();
1962  TIBL3RX.clear();
1963  TIBL4RX.clear();
1964  TIBL1RY.clear();
1965  TIBL2RY.clear();
1966  TIBL3RY.clear();
1967  TIBL4RY.clear();
1968  TIBL1SX.clear();
1969  TIBL2SX.clear();
1970  TIBL3SX.clear();
1971  TIBL4SX.clear();
1972  TIBL1SY.clear();
1973  TIBL2SY.clear();
1974  TIBL3SY.clear();
1975  TIBL4SY.clear();
1976 
1977  TOBL1RX.clear();
1978  TOBL2RX.clear();
1979  TOBL3RX.clear();
1980  TOBL4RX.clear();
1981  TOBL1RY.clear();
1982  TOBL2RY.clear();
1983  TOBL3RY.clear();
1984  TOBL4RY.clear();
1985  TOBL1SX.clear();
1986  TOBL2SX.clear();
1987  TOBL3SX.clear();
1988  TOBL4SX.clear();
1989  TOBL1SY.clear();
1990  TOBL2SY.clear();
1991  TOBL3SY.clear();
1992  TOBL4SY.clear();
1993 
1994  TIDW1RX.clear();
1995  TIDW2RX.clear();
1996  TIDW3RX.clear();
1997  TIDW1RY.clear();
1998  TIDW2RY.clear();
1999  TIDW3RY.clear();
2000  TIDW1SX.clear();
2001  TIDW2SX.clear();
2002  TIDW3SX.clear();
2003  TIDW1SY.clear();
2004  TIDW2SY.clear();
2005  TIDW3SY.clear();
2006 
2007  TECW1RX.clear();
2008  TECW2RX.clear();
2009  TECW3RX.clear();
2010  TECW4RX.clear();
2011  TECW5RX.clear();
2012  TECW6RX.clear();
2013  TECW7RX.clear();
2014  TECW8RX.clear();
2015  TECW1RY.clear();
2016  TECW2RY.clear();
2017  TECW3RY.clear();
2018  TECW4RY.clear();
2019  TECW5RY.clear();
2020  TECW6RY.clear();
2021  TECW7RY.clear();
2022  TECW8RY.clear();
2023  TECW1SX.clear();
2024  TECW2SX.clear();
2025  TECW3SX.clear();
2026  TECW4SX.clear();
2027  TECW5SX.clear();
2028  TECW6SX.clear();
2029  TECW7SX.clear();
2030  TECW8SX.clear();
2031  TECW1SY.clear();
2032  TECW2SY.clear();
2033  TECW3SY.clear();
2034  TECW4SY.clear();
2035  TECW5SY.clear();
2036  TECW6SY.clear();
2037  TECW7SY.clear();
2038  TECW8SY.clear();
2039 
2040  BRL1RX.clear();
2041  BRL1RY.clear();
2042  BRL1SX.clear();
2043  BRL1SY.clear();
2044  BRL2RX.clear();
2045  BRL2RY.clear();
2046  BRL2SX.clear();
2047  BRL2SY.clear();
2048  BRL3RX.clear();
2049  BRL3RY.clear();
2050  BRL3SX.clear();
2051  BRL3SY.clear();
2052 
2053  FWD1pRX.clear();
2054  FWD1pRY.clear();
2055  FWD1pSX.clear();
2056  FWD1pSY.clear();
2057  FWD1nRX.clear();
2058  FWD1nRY.clear();
2059  FWD1nSX.clear();
2060  FWD1nSY.clear();
2061  FWD2pRX.clear();
2062  FWD2pRY.clear();
2063  FWD2pSX.clear();
2064  FWD2pSY.clear();
2065  FWD2nRX.clear();
2066  FWD2nRY.clear();
2067  FWD2nSX.clear();
2068  FWD2nSY.clear();
2069 
2070  //muon clear
2071  DTRHD.clear();
2072  DTSHD.clear();
2073 
2074  CSCRHPHI.clear();
2075  CSCRHPERP.clear();
2076  CSCSHPHI.clear();
2077 
2078  RPCRHX.clear();
2079  RPCSHX.clear();
2080 
2081  return;
2082 }
2083 
2084 //needed by to do the residual for matched hits in SiStrip
2085 std::pair<LocalPoint,LocalVector>
2087  const StripGeomDetUnit* stripDet,
2088  const BoundPlane& plane)
2089 {
2090 
2091  const StripTopology& topol = stripDet->specificTopology();
2092  GlobalPoint globalpos= stripDet->surface().toGlobal(hit.localPosition());
2093  LocalPoint localHit = plane.toLocal(globalpos);
2094  //track direction
2095  LocalVector locdir=hit.localDirection();
2096  //rotate track in new frame
2097 
2098  GlobalVector globaldir= stripDet->surface().toGlobal(locdir);
2099  LocalVector dir=plane.toLocal(globaldir);
2100  float scale = -localHit.z() / dir.z();
2101 
2102  LocalPoint projectedPos = localHit + scale*dir;
2103 
2104  float selfAngle = topol.stripAngle( topol.strip( hit.localPosition()));
2105 
2106  // vector along strip in hit frame
2107  LocalVector stripDir( sin(selfAngle), cos(selfAngle), 0);
2108 
2109  LocalVector
2110  localStripDir(plane.toLocal(stripDet->surface().toGlobal(stripDir)));
2111 
2112  return std::pair<LocalPoint,LocalVector>( projectedPos, localStripDir);
2113 }
2114 
2115 // Return a map between DTRecHit1DPair and wireId
2116 std::map<DTWireId, std::vector<DTRecHit1DPair> >
2118  dt1DRecHitPairs) {
2119  std::map<DTWireId, std::vector<DTRecHit1DPair> > ret;
2120 
2121  for(DTRecHitCollection::const_iterator rechit = dt1DRecHitPairs->begin();
2122  rechit != dt1DRecHitPairs->end(); rechit++) {
2123  ret[(*rechit).wireId()].push_back(*rechit);
2124  }
2125 
2126  return ret;
2127 }
2128 
2129 // Compute SimHit distance from wire (cm)
2131  DTWireId wireId,
2132  const PSimHit& hit) {
2133  float xwire = layer->specificTopology().wirePosition(wireId.wire());
2134  LocalPoint entryP = hit.entryPoint();
2135  LocalPoint exitP = hit.exitPoint();
2136  float xEntry = entryP.x()-xwire;
2137  float xExit = exitP.x()-xwire;
2138 
2139  //FIXME: check...
2140  return fabs(xEntry - (entryP.z()*(xExit-xEntry))/(exitP.z()-entryP.z()));
2141 }
2142 
2143 // Find the RecHit closest to the muon SimHit
2144 template <typename type>
2145 const type*
2147  DTWireId wireId,
2148  const std::vector<type>& recHits,
2149  const float simHitDist) {
2150  float res = 99999;
2151  const type* theBestRecHit = 0;
2152  // Loop over RecHits within the cell
2153  for(typename std::vector<type>::const_iterator recHit = recHits.begin();
2154  recHit != recHits.end();
2155  recHit++) {
2156  float distTmp = recHitDistFromWire(*recHit, layer);
2157  if(fabs(distTmp-simHitDist) < res) {
2158  res = fabs(distTmp-simHitDist);
2159  theBestRecHit = &(*recHit);
2160  }
2161  } // End of loop over RecHits within the cell
2162 
2163  return theBestRecHit;
2164 }
2165 
2166 // Compute the distance from wire (cm) of a hits in a DTRecHit1DPair
2167 float
2169  const DTLayer* layer) {
2170  // Compute the rechit distance from wire
2171  return fabs(hitPair.localPosition(DTEnums::Left).x() -
2172  hitPair.localPosition(DTEnums::Right).x())/2.;
2173 }
2174 
2175 // Compute the distance from wire (cm) of a hits in a DTRecHit1D
2176 float
2178  const DTLayer* layer) {
2179  return fabs(recHit.localPosition().x() -
2180  layer->specificTopology().wirePosition(recHit.wireId().wire()));
2181 }
2182 
2183 template <typename type>
2185  std::map<DTWireId, std::vector<PSimHit> >
2186  simHitsPerWire,
2187  std::map<DTWireId, std::vector<type> >
2188  recHitsPerWire,
2189  int step) {
2190 
2191  int nDt = 0;
2192  // Loop over cells with a muon SimHit
2193  for(std::map<DTWireId, std::vector<PSimHit> >::const_iterator wireAndSHits =
2194  simHitsPerWire.begin();
2195  wireAndSHits != simHitsPerWire.end();
2196  wireAndSHits++) {
2197  DTWireId wireId = (*wireAndSHits).first;
2198  std::vector<PSimHit> simHitsInCell = (*wireAndSHits).second;
2199 
2200  // Get the layer
2201  const DTLayer* layer = dtGeom->layer(wireId);
2202 
2203  // Look for a mu hit in the cell
2204  const PSimHit* muSimHit = DTHitQualityUtils::findMuSimHit(simHitsInCell);
2205  if (muSimHit==0) {
2206  continue; // Skip this cell
2207  }
2208 
2209  // Find the distance of the simhit from the wire
2210  float simHitWireDist = simHitDistFromWire(layer, wireId, *muSimHit);
2211  // Skip simhits out of the cell
2212  if(simHitWireDist>2.1) {
2213  continue; // Skip this cell
2214  }
2215  //GlobalPoint simHitGlobalPos = layer->toGlobal(muSimHit->localPosition());
2216 
2217  // Look for RecHits in the same cell
2218  if(recHitsPerWire.find(wireId) == recHitsPerWire.end()) {
2219  continue; // No RecHit found in this cell
2220  } else {
2221 
2222  // vector<type> recHits = (*wireAndRecHits).second;
2223  std::vector<type> recHits = recHitsPerWire[wireId];
2224 
2225  // Find the best RecHit
2226  const type* theBestRecHit =
2227  findBestRecHit(layer, wireId, recHits, simHitWireDist);
2228 
2229  float recHitWireDist = recHitDistFromWire(*theBestRecHit, layer);
2230 
2231  ++nDt;
2232 
2233  DTRHD.push_back(recHitWireDist);
2234  DTSHD.push_back(simHitWireDist);
2235 
2236  } // find rechits
2237  } // loop over simhits
2238 
2239  return nDt;
2240 }
2241 
2242 void
2244  const CSCRecHit2D & recHit,
2245  const CSCLayer * layer,
2246  int chamberType) {
2247  GlobalPoint simHitPos = layer->toGlobal(simHit.localPosition());
2248  GlobalPoint recHitPos = layer->toGlobal(recHit.localPosition());
2249 
2250  CSCRHPHI.push_back(recHitPos.phi());
2251  CSCRHPERP.push_back(recHitPos.perp());
2252  CSCSHPHI.push_back(simHitPos.phi());
2253 }
2254 
2255 //define this as a plug-in
2256 //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:407
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:73
void getAllProvenance(std::vector< Provenance const * > &provenances) const
Definition: Event.cc:71
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:66
unsigned int layer() const
layer id
Definition: TOBDetId.h:39
const std::string & label
Definition: MVAComputer.cc:186
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:49
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:63
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:57
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 &)
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 &)
virtual LocalPoint localPosition() const
unsigned int layer() const
layer id
Definition: PXBDetId.h:35
std::map< int, edm::PSimHitContainer > theMap
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:46
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:84
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:28
LocalPoint toLocal(const GlobalPoint &gp) const
T z() const
Definition: PV3DBase.h:58
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:355
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 int verbosity
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:70
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
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)
DetId geographicalId() const
dbl *** dir
Definition: mlp_gen.cc:35
static std::map< DTWireId, edm::PSimHitContainer > mapSimHitsPerWire(const edm::PSimHitContainer &simhits)
LocalPoint localPosition() const
Definition: CSCRecHit2D.h:45
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 &)
const GlobalPoint & getPosition() const
bool isForward() const
Definition: RPCRoll.cc:98
T x() const
Definition: PV3DBase.h:56
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
virtual const BoundPlane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
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