GlobalRecHits/src/GlobalRecHitsProducer.cc

0001 /** \file GlobalRecHitsProducer.cc
0002  *
0003  *  See header file for description of class
0004  *
0005  *  \author M. Strang SUNY-Buffalo
0006  */
0007
0008 #include "Validation/GlobalRecHits/interface/GlobalRecHitsProducer.h"
0009 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0010 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0011 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0012 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0013
0014 GlobalRecHitsProducer::GlobalRecHitsProducer(const edm::ParameterSet& iPSet)
0015     : fName(""),
0016       verbosity(0),
0017       frequency(0),
0018       label(""),
0019       getAllProvenances(false),
0020       printProvenanceInfo(false),
0021       trackerHitAssociatorConfig_(iPSet, consumesCollector()),
0022       caloGeomToken_(esConsumes()),
0023       tTopoToken_(esConsumes()),
0024       tGeomToken_(esConsumes()),
0025       dtGeomToken_(esConsumes()),
0026       cscGeomToken_(esConsumes()),
0027       rpcGeomToken_(esConsumes()),
0028       count(0) {
0029   std::string MsgLoggerCat = "GlobalRecHitsProducer_GlobalRecHitsProducer";
0030
0031   // get information from parameter set
0032   fName = iPSet.getUntrackedParameter<std::string>("Name");
0033   verbosity = iPSet.getUntrackedParameter<int>("Verbosity");
0034   frequency = iPSet.getUntrackedParameter<int>("Frequency");
0035   label = iPSet.getParameter<std::string>("Label");
0036   edm::ParameterSet m_Prov = iPSet.getParameter<edm::ParameterSet>("ProvenanceLookup");
0037   getAllProvenances = m_Prov.getUntrackedParameter<bool>("GetAllProvenances");
0038   printProvenanceInfo = m_Prov.getUntrackedParameter<bool>("PrintProvenanceInfo");
0039
0040   //get Labels to use to extract information
0041   ECalEBSrc_ = iPSet.getParameter<edm::InputTag>("ECalEBSrc");
0042   ECalUncalEBSrc_ = iPSet.getParameter<edm::InputTag>("ECalUncalEBSrc");
0043   ECalEESrc_ = iPSet.getParameter<edm::InputTag>("ECalEESrc");
0044   ECalUncalEESrc_ = iPSet.getParameter<edm::InputTag>("ECalUncalEESrc");
0045   ECalESSrc_ = iPSet.getParameter<edm::InputTag>("ECalESSrc");
0046   HCalSrc_ = iPSet.getParameter<edm::InputTag>("HCalSrc");
0047   SiStripSrc_ = iPSet.getParameter<edm::InputTag>("SiStripSrc");
0048   SiPxlSrc_ = iPSet.getParameter<edm::InputTag>("SiPxlSrc");
0049   MuDTSrc_ = iPSet.getParameter<edm::InputTag>("MuDTSrc");
0050   MuDTSimSrc_ = iPSet.getParameter<edm::InputTag>("MuDTSimSrc");
0051   MuCSCSrc_ = iPSet.getParameter<edm::InputTag>("MuCSCSrc");
0052   MuRPCSrc_ = iPSet.getParameter<edm::InputTag>("MuRPCSrc");
0053   MuRPCSimSrc_ = iPSet.getParameter<edm::InputTag>("MuRPCSimSrc");
0054
0055   // fix for consumes
0056   ECalUncalEBSrc_Token_ = consumes<EBUncalibratedRecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalUncalEBSrc"));
0057   ECalUncalEESrc_Token_ = consumes<EEUncalibratedRecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalUncalEESrc"));
0058   ECalEBSrc_Token_ = consumes<EBRecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalEBSrc"));
0059   ECalEESrc_Token_ = consumes<EERecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalEESrc"));
0060   ECalESSrc_Token_ = consumes<ESRecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalESSrc"));
0061   HCalSrc_Token_ = consumes<edm::PCaloHitContainer>(iPSet.getParameter<edm::InputTag>("HCalSrc"));
0062   SiStripSrc_Token_ = consumes<SiStripMatchedRecHit2DCollection>(iPSet.getParameter<edm::InputTag>("SiStripSrc"));
0063   SiPxlSrc_Token_ = consumes<SiPixelRecHitCollection>(iPSet.getParameter<edm::InputTag>("SiPxlSrc"));
0064
0065   MuDTSrc_Token_ = consumes<DTRecHitCollection>(iPSet.getParameter<edm::InputTag>("MuDTSrc"));
0066   MuDTSimSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("MuDTSimSrc"));
0067
0068   MuCSCSrc_Token_ = consumes<CSCRecHit2DCollection>(iPSet.getParameter<edm::InputTag>("MuCSCSrc"));
0069   MuCSCHits_Token_ = consumes<CrossingFrame<PSimHit>>(
0070       edm::InputTag(std::string("mix"), iPSet.getParameter<std::string>("hitsProducer") + std::string("MuonCSCHits")));
0071
0072   MuRPCSrc_Token_ = consumes<RPCRecHitCollection>(iPSet.getParameter<edm::InputTag>("MuRPCSrc"));
0073   MuRPCSimSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("MuRPCSimSrc"));
0074
0075   EBHits_Token_ = consumes<CrossingFrame<PCaloHit>>(
0076       edm::InputTag(std::string("mix"), iPSet.getParameter<std::string>("hitsProducer") + std::string("EcalHitsEB")));
0077   EEHits_Token_ =
0078       consumes<CrossingFrame<PCaloHit>>(edm::InputTag(std::string("mix"), iPSet.getParameter<std::string>("hitsProduc\
0079 er") + std::string("EcalHitsEE")));
0080
0081   // use value of first digit to determine default output level (inclusive)
0082   // 0 is none, 1 is basic, 2 is fill output, 3 is gather output
0083   verbosity %= 10;
0084
0085   // create persistent object
0086   produces<PGlobalRecHit>(label);
0087
0088   // print out Parameter Set information being used
0089   if (verbosity >= 0) {
0090     edm::LogInfo(MsgLoggerCat) << "\n===============================\n"
0091                                << "Initialized as EDProducer with parameter values:\n"
0092                                << "    Name           = " << fName << "\n"
0093                                << "    Verbosity      = " << verbosity << "\n"
0094                                << "    Frequency      = " << frequency << "\n"
0095                                << "    Label          = " << label << "\n"
0096                                << "    GetProv        = " << getAllProvenances << "\n"
0097                                << "    PrintProv      = " << printProvenanceInfo << "\n"
0098                                << "    ECalEBSrc      = " << ECalEBSrc_.label() << ":" << ECalEBSrc_.instance() << "\n"
0099                                << "    ECalUncalEBSrc = " << ECalUncalEBSrc_.label() << ":"
0100                                << ECalUncalEBSrc_.instance() << "\n"
0101                                << "    ECalEESrc      = " << ECalEESrc_.label() << ":" << ECalUncalEESrc_.instance()
0102                                << "\n"
0103                                << "    ECalUncalEESrc = " << ECalUncalEESrc_.label() << ":" << ECalEESrc_.instance()
0104                                << "\n"
0105                                << "    ECalESSrc      = " << ECalESSrc_.label() << ":" << ECalESSrc_.instance() << "\n"
0106                                << "    HCalSrc        = " << HCalSrc_.label() << ":" << HCalSrc_.instance() << "\n"
0107                                << "    SiStripSrc     = " << SiStripSrc_.label() << ":" << SiStripSrc_.instance()
0108                                << "\n"
0109                                << "    SiPixelSrc     = " << SiPxlSrc_.label() << ":" << SiPxlSrc_.instance() << "\n"
0110                                << "    MuDTSrc        = " << MuDTSrc_.label() << ":" << MuDTSrc_.instance() << "\n"
0111                                << "    MuDTSimSrc     = " << MuDTSimSrc_.label() << ":" << MuDTSimSrc_.instance()
0112                                << "\n"
0113                                << "    MuCSCSrc       = " << MuCSCSrc_.label() << ":" << MuCSCSrc_.instance() << "\n"
0114                                << "    MuRPCSrc       = " << MuRPCSrc_.label() << ":" << MuRPCSrc_.instance() << "\n"
0115                                << "    MuRPCSimSrc    = " << MuRPCSimSrc_.label() << ":" << MuRPCSimSrc_.instance()
0116                                << "\n"
0117                                << "===============================\n";
0118   }
0119 }
0120
0121 GlobalRecHitsProducer::~GlobalRecHitsProducer() {}
0122
0123 void GlobalRecHitsProducer::beginJob() {
0124   std::string MsgLoggerCat = "GlobalRecHitsProducer_beginJob";
0125
0126   // clear storage vectors
0127   clear();
0128   return;
0129 }
0130
0131 void GlobalRecHitsProducer::endJob() {
0132   std::string MsgLoggerCat = "GlobalRecHitsProducer_endJob";
0133   if (verbosity >= 0)
0134     edm::LogInfo(MsgLoggerCat) << "Terminating having processed " << count << " events.";
0135   return;
0136 }
0137
0138 void GlobalRecHitsProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0139   std::string MsgLoggerCat = "GlobalRecHitsProducer_produce";
0140
0141   // keep track of number of events processed
0142   ++count;
0143
0144   // get event id information
0145   edm::RunNumber_t nrun = iEvent.id().run();
0146   edm::EventNumber_t nevt = iEvent.id().event();
0147
0148   if (verbosity > 0) {
0149     edm::LogInfo(MsgLoggerCat) << "Processing run " << nrun << ", event " << nevt << " (" << count << " events total)";
0150   } else if (verbosity == 0) {
0151     if (nevt % frequency == 0 || nevt == 1) {
0152       edm::LogInfo(MsgLoggerCat) << "Processing run " << nrun << ", event " << nevt << " (" << count
0153                                  << " events total)";
0154     }
0155   }
0156
0157   // clear event holders
0158   clear();
0159
0160   // look at information available in the event
0161   if (getAllProvenances) {
0162     std::vector<const edm::StableProvenance*> AllProv;
0163     iEvent.getAllStableProvenance(AllProv);
0164
0165     if (verbosity >= 0)
0166       edm::LogInfo(MsgLoggerCat) << "Number of Provenances = " << AllProv.size();
0167
0168     if (printProvenanceInfo && (verbosity >= 0)) {
0169       TString eventout("\nProvenance info:\n");
0170
0171       for (unsigned int i = 0; i < AllProv.size(); ++i) {
0172         eventout += "\n       ******************************";
0173         eventout += "\n       Module       : ";
0174         //eventout += (AllProv[i]->product).moduleLabel();
0175         eventout += AllProv[i]->moduleLabel();
0176         eventout += "\n       ProductID    : ";
0177         //eventout += (AllProv[i]->product).productID_.id_;
0178         eventout += AllProv[i]->productID().id();
0179         eventout += "\n       ClassName    : ";
0180         //eventout += (AllProv[i]->product).fullClassName_;
0181         eventout += AllProv[i]->className();
0182         eventout += "\n       InstanceName : ";
0183         //eventout += (AllProv[i]->product).productInstanceName_;
0184         eventout += AllProv[i]->productInstanceName();
0185         eventout += "\n       BranchName   : ";
0186         //eventout += (AllProv[i]->product).branchName_;
0187         eventout += AllProv[i]->branchName();
0188       }
0189       eventout += "\n       ******************************\n";
0190       edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0191       printProvenanceInfo = false;
0192     }
0193     getAllProvenances = false;
0194   }
0195
0196   // call fill functions
0197   // gather Ecal information from event
0198   fillECal(iEvent, iSetup);
0199   // gather Hcal information from event
0200   fillHCal(iEvent, iSetup);
0201   // gather Track information from event
0202   fillTrk(iEvent, iSetup);
0203   // gather Muon information from event
0204   fillMuon(iEvent, iSetup);
0205
0206   if (verbosity > 0)
0207     edm::LogInfo(MsgLoggerCat) << "Done gathering data from event.";
0208
0209   // produce object to put into event
0210   std::unique_ptr<PGlobalRecHit> pOut(new PGlobalRecHit);
0211
0212   if (verbosity > 2)
0213     edm::LogInfo(MsgLoggerCat) << "Saving event contents:";
0214
0215   // call store functions
0216   // store ECal information in produce
0217   storeECal(*pOut);
0218   // store HCal information in produce
0219   storeHCal(*pOut);
0220   // store Track information in produce
0221   storeTrk(*pOut);
0222   // store Muon information in produce
0223   storeMuon(*pOut);
0224
0225   // store information in event
0226   iEvent.put(std::move(pOut), label);
0227
0228   return;
0229 }
0230
0231 void GlobalRecHitsProducer::fillECal(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0232   std::string MsgLoggerCat = "GlobalRecHitsProducer_fillECal";
0233
0234   TString eventout;
0235   if (verbosity > 0)
0236     eventout = "\nGathering info:";
0237
0238   // extract crossing frame from event
0239   //edm::Handle<CrossingFrame> crossingFrame;
0240   edm::Handle<CrossingFrame<PCaloHit>> crossingFrame;
0241   //iEvent.getByType(crossingFrame);
0242   //if (!crossingFrame.isValid()) {
0243   //  edm::LogWarning(MsgLoggerCat)
0244   //    << "Unable to find crossingFrame in event!";
0245   //  return;
0246   //}
0247
0248   ////////////////////////
0249   //extract EB information
0250   ////////////////////////
0251   edm::Handle<EBUncalibratedRecHitCollection> EcalUncalibRecHitEB;
0252   iEvent.getByToken(ECalUncalEBSrc_Token_, EcalUncalibRecHitEB);
0253   if (!EcalUncalibRecHitEB.isValid()) {
0254     edm::LogWarning(MsgLoggerCat) << "Unable to find EcalUncalRecHitEB in event!";
0255     return;
0256   }
0257
0258   edm::Handle<EBRecHitCollection> EcalRecHitEB;
0259   iEvent.getByToken(ECalEBSrc_Token_, EcalRecHitEB);
0260   if (!EcalRecHitEB.isValid()) {
0261     edm::LogWarning(MsgLoggerCat) << "Unable to find EcalRecHitEB in event!";
0262     return;
0263   }
0264
0265   // loop over simhits
0266   iEvent.getByToken(EBHits_Token_, crossingFrame);
0267   if (!crossingFrame.isValid()) {
0268     edm::LogWarning(MsgLoggerCat) << "Unable to find cal barrel crossingFrame in event!";
0269     return;
0270   }
0271   //std::unique_ptr<MixCollection<PCaloHit> >
0272   //  barrelHits(new MixCollection<PCaloHit>
0273   //           (crossingFrame.product(), barrelHitsName));
0274   std::unique_ptr<MixCollection<PCaloHit>> barrelHits(new MixCollection<PCaloHit>(crossingFrame.product()));
0275
0276   // keep track of sum of simhit energy in each crystal
0277   MapType ebSimMap;
0278   for (MixCollection<PCaloHit>::MixItr hitItr = barrelHits->begin(); hitItr != barrelHits->end(); ++hitItr) {
0279     EBDetId ebid = EBDetId(hitItr->id());
0280
0281     uint32_t crystid = ebid.rawId();
0282     ebSimMap[crystid] += hitItr->energy();
0283   }
0284
0285   int nEBRecHits = 0;
0286   // loop over RecHits
0287   const EBUncalibratedRecHitCollection* EBUncalibRecHit = EcalUncalibRecHitEB.product();
0288   const EBRecHitCollection* EBRecHit = EcalRecHitEB.product();
0289
0290   for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EBUncalibRecHit->begin();
0291        uncalibRecHit != EBUncalibRecHit->end();
0292        ++uncalibRecHit) {
0293     EBDetId EBid = EBDetId(uncalibRecHit->id());
0294
0295     EcalRecHitCollection::const_iterator myRecHit = EBRecHit->find(EBid);
0296
0297     if (myRecHit != EBRecHit->end()) {
0298       ++nEBRecHits;
0299       EBRE.push_back(myRecHit->energy());
0300       EBSHE.push_back(ebSimMap[EBid.rawId()]);
0301     }
0302   }
0303
0304   if (verbosity > 1) {
0305     eventout += "\n          Number of EBRecHits collected:............ ";
0306     eventout += nEBRecHits;
0307   }
0308
0309   ////////////////////////
0310   //extract EE information
0311   ////////////////////////
0312   edm::Handle<EEUncalibratedRecHitCollection> EcalUncalibRecHitEE;
0313   iEvent.getByToken(ECalUncalEESrc_Token_, EcalUncalibRecHitEE);
0314   if (!EcalUncalibRecHitEE.isValid()) {
0315     edm::LogWarning(MsgLoggerCat) << "Unable to find EcalUncalRecHitEE in event!";
0316     return;
0317   }
0318
0319   edm::Handle<EERecHitCollection> EcalRecHitEE;
0320   iEvent.getByToken(ECalEESrc_Token_, EcalRecHitEE);
0321   if (!EcalRecHitEE.isValid()) {
0322     edm::LogWarning(MsgLoggerCat) << "Unable to find EcalRecHitEE in event!";
0323     return;
0324   }
0325
0326   // loop over simhits
0327   iEvent.getByToken(EEHits_Token_, crossingFrame);
0328   if (!crossingFrame.isValid()) {
0329     edm::LogWarning(MsgLoggerCat) << "Unable to find cal endcap crossingFrame in event!";
0330     return;
0331   }
0332   //std::unique_ptr<MixCollection<PCaloHit> >
0333   //  endcapHits(new MixCollection<PCaloHit>
0334   //           (crossingFrame.product(), endcapHitsName));
0335   std::unique_ptr<MixCollection<PCaloHit>> endcapHits(new MixCollection<PCaloHit>(crossingFrame.product()));
0336
0337   // keep track of sum of simhit energy in each crystal
0338   MapType eeSimMap;
0339   for (MixCollection<PCaloHit>::MixItr hitItr = endcapHits->begin(); hitItr != endcapHits->end(); ++hitItr) {
0340     EEDetId eeid = EEDetId(hitItr->id());
0341
0342     uint32_t crystid = eeid.rawId();
0343     eeSimMap[crystid] += hitItr->energy();
0344   }
0345
0346   int nEERecHits = 0;
0347   // loop over RecHits
0348   const EEUncalibratedRecHitCollection* EEUncalibRecHit = EcalUncalibRecHitEE.product();
0349   const EERecHitCollection* EERecHit = EcalRecHitEE.product();
0350
0351   for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EEUncalibRecHit->begin();
0352        uncalibRecHit != EEUncalibRecHit->end();
0353        ++uncalibRecHit) {
0354     EEDetId EEid = EEDetId(uncalibRecHit->id());
0355
0356     EcalRecHitCollection::const_iterator myRecHit = EERecHit->find(EEid);
0357
0358     if (myRecHit != EERecHit->end()) {
0359       ++nEERecHits;
0360       EERE.push_back(myRecHit->energy());
0361       EESHE.push_back(eeSimMap[EEid.rawId()]);
0362     }
0363   }
0364
0365   if (verbosity > 1) {
0366     eventout += "\n          Number of EERecHits collected:............ ";
0367     eventout += nEERecHits;
0368   }
0369
0370   ////////////////////////
0371   //extract ES information
0372   ////////////////////////
0373   edm::Handle<ESRecHitCollection> EcalRecHitES;
0374   iEvent.getByToken(ECalESSrc_Token_, EcalRecHitES);
0375   if (!EcalRecHitES.isValid()) {
0376     edm::LogWarning(MsgLoggerCat) << "Unable to find EcalRecHitES in event!";
0377     return;
0378   }
0379
0380   // loop over simhits
0381   iEvent.getByToken(ESHits_Token_, crossingFrame);
0382   if (!crossingFrame.isValid()) {
0383     edm::LogWarning(MsgLoggerCat) << "Unable to find cal preshower crossingFrame in event!";
0384     return;
0385   }
0386   //std::unique_ptr<MixCollection<PCaloHit> >
0387   //  preshowerHits(new MixCollection<PCaloHit>
0388   //           (crossingFrame.product(), preshowerHitsName));
0389   std::unique_ptr<MixCollection<PCaloHit>> preshowerHits(new MixCollection<PCaloHit>(crossingFrame.product()));
0390
0391   // keep track of sum of simhit energy in each crystal
0392   MapType esSimMap;
0393   for (MixCollection<PCaloHit>::MixItr hitItr = preshowerHits->begin(); hitItr != preshowerHits->end(); ++hitItr) {
0394     ESDetId esid = ESDetId(hitItr->id());
0395
0396     uint32_t crystid = esid.rawId();
0397     esSimMap[crystid] += hitItr->energy();
0398   }
0399
0400   int nESRecHits = 0;
0401   // loop over RecHits
0402   const ESRecHitCollection* ESRecHit = EcalRecHitES.product();
0403   for (EcalRecHitCollection::const_iterator recHit = ESRecHit->begin(); recHit != ESRecHit->end(); ++recHit) {
0404     ESDetId ESid = ESDetId(recHit->id());
0405
0406     ++nESRecHits;
0407     ESRE.push_back(recHit->energy());
0408     ESSHE.push_back(esSimMap[ESid.rawId()]);
0409   }
0410
0411   if (verbosity > 1) {
0412     eventout += "\n          Number of ESRecHits collected:............ ";
0413     eventout += nESRecHits;
0414   }
0415
0416   if (verbosity > 0)
0417     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0418
0419   return;
0420 }
0421
0422 void GlobalRecHitsProducer::storeECal(PGlobalRecHit& product) {
0423   std::string MsgLoggerCat = "GlobalRecHitsProducer_storeECal";
0424
0425   if (verbosity > 2) {
0426     TString eventout("\n         nEBRecHits     = ");
0427     eventout += EBRE.size();
0428     for (unsigned int i = 0; i < EBRE.size(); ++i) {
0429       eventout += "\n      (RE, SHE) = (";
0430       eventout += EBRE[i];
0431       eventout += ", ";
0432       eventout += EBSHE[i];
0433       eventout += ")";
0434     }
0435     eventout += "\n         nEERecHits     = ";
0436     eventout += EERE.size();
0437     for (unsigned int i = 0; i < EERE.size(); ++i) {
0438       eventout += "\n      (RE, SHE) = (";
0439       eventout += EERE[i];
0440       eventout += ", ";
0441       eventout += EESHE[i];
0442       eventout += ")";
0443     }
0444     eventout += "\n         nESRecHits     = ";
0445     eventout += ESRE.size();
0446     for (unsigned int i = 0; i < ESRE.size(); ++i) {
0447       eventout += "\n      (RE, SHE) = (";
0448       eventout += ESRE[i];
0449       eventout += ", ";
0450       eventout += ESSHE[i];
0451       eventout += ")";
0452     }
0453     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0454   }
0455
0456   product.putEBCalRecHits(EBRE, EBSHE);
0457   product.putEECalRecHits(EERE, EESHE);
0458   product.putESCalRecHits(ESRE, ESSHE);
0459
0460   return;
0461 }
0462
0463 void GlobalRecHitsProducer::fillHCal(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0464   std::string MsgLoggerCat = "GlobalRecHitsProducer_fillHCal";
0465
0466   TString eventout;
0467   if (verbosity > 0)
0468     eventout = "\nGathering info:";
0469
0470   // get geometry
0471   const auto& geometry = iSetup.getHandle(caloGeomToken_);
0472   if (!geometry.isValid()) {
0473     edm::LogWarning(MsgLoggerCat) << "Unable to find CaloGeometry in event!";
0474     return;
0475   }
0476
0477   ///////////////////////
0478   // extract simhit info
0479   //////////////////////
0480   edm::Handle<edm::PCaloHitContainer> hcalHits;
0481   iEvent.getByToken(HCalSrc_Token_, hcalHits);
0482   if (!hcalHits.isValid()) {
0483     edm::LogWarning(MsgLoggerCat) << "Unable to find hcalHits in event!";
0484     return;
0485   }
0486   const edm::PCaloHitContainer* simhitResult = hcalHits.product();
0487
0488   MapType fHBEnergySimHits;
0489   MapType fHEEnergySimHits;
0490   MapType fHOEnergySimHits;
0491   MapType fHFEnergySimHits;
0492   for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin(); simhits != simhitResult->end();
0493        ++simhits) {
0494     HcalDetId detId(simhits->id());
0495     uint32_t cellid = detId.rawId();
0496
0497     if (detId.subdet() == sdHcalBrl) {
0498       fHBEnergySimHits[cellid] += simhits->energy();
0499     }
0500     if (detId.subdet() == sdHcalEC) {
0501       fHEEnergySimHits[cellid] += simhits->energy();
0502     }
0503     if (detId.subdet() == sdHcalOut) {
0504       fHOEnergySimHits[cellid] += simhits->energy();
0505     }
0506     if (detId.subdet() == sdHcalFwd) {
0507       fHFEnergySimHits[cellid] += simhits->energy();
0508     }
0509   }
0510
0511   // max values to be used (HO is found in HB)
0512   Double_t maxHBEnergy = 0.;
0513   Double_t maxHEEnergy = 0.;
0514   Double_t maxHFEnergy = 0.;
0515
0516   Double_t maxHBPhi = -1000.;
0517   Double_t maxHEPhi = -1000.;
0518   Double_t maxHOPhi = -1000.;
0519   Double_t maxHFPhi = -1000.;
0520
0521   Double_t maxHBEta = -1000.;
0522   Double_t maxHEEta = -1000.;
0523   Double_t maxHOEta = -1000.;
0524   Double_t maxHFEta = -1000.;
0525
0526   Double_t PI = 3.141592653589;
0527
0528   ////////////////////////
0529   // get HBHE information
0530   ///////////////////////
0531   std::vector<edm::Handle<HBHERecHitCollection>> hbhe;
0532   iEvent.getManyByType(hbhe);
0533   if (!hbhe[0].isValid()) {
0534     edm::LogWarning(MsgLoggerCat) << "Unable to find any HBHERecHitCollections in event!";
0535     return;
0536   }
0537   std::vector<edm::Handle<HBHERecHitCollection>>::iterator ihbhe;
0538   const CaloGeometry* geo = geometry.product();
0539
0540   int iHB = 0;
0541   int iHE = 0;
0542   for (ihbhe = hbhe.begin(); ihbhe != hbhe.end(); ++ihbhe) {
0543     // find max values
0544     for (HBHERecHitCollection::const_iterator jhbhe = (*ihbhe)->begin(); jhbhe != (*ihbhe)->end(); ++jhbhe) {
0545       HcalDetId cell(jhbhe->id());
0546
0547       if (cell.subdet() == sdHcalBrl) {
0548         const HcalGeometry* cellGeometry =
0549             dynamic_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, cell.subdet()));
0550         double fEta = cellGeometry->getPosition(cell).eta();
0551         double fPhi = cellGeometry->getPosition(cell).phi();
0552         if ((jhbhe->energy()) > maxHBEnergy) {
0553           maxHBEnergy = jhbhe->energy();
0554           maxHBPhi = fPhi;
0555           maxHOPhi = maxHBPhi;
0556           maxHBEta = fEta;
0557           maxHOEta = maxHBEta;
0558         }
0559       }
0560
0561       if (cell.subdet() == sdHcalEC) {
0562         const HcalGeometry* cellGeometry =
0563             dynamic_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, cell.subdet()));
0564         double fEta = cellGeometry->getPosition(cell).eta();
0565         double fPhi = cellGeometry->getPosition(cell).phi();
0566         if ((jhbhe->energy()) > maxHEEnergy) {
0567           maxHEEnergy = jhbhe->energy();
0568           maxHEPhi = fPhi;
0569           maxHEEta = fEta;
0570         }
0571       }
0572     }  // end find max values
0573
0574     for (HBHERecHitCollection::const_iterator jhbhe = (*ihbhe)->begin(); jhbhe != (*ihbhe)->end(); ++jhbhe) {
0575       HcalDetId cell(jhbhe->id());
0576
0577       if (cell.subdet() == sdHcalBrl) {
0578         ++iHB;
0579
0580         const HcalGeometry* cellGeometry =
0581             dynamic_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, cell.subdet()));
0582         double fEta = cellGeometry->getPosition(cell).eta();
0583         double fPhi = cellGeometry->getPosition(cell).phi();
0584
0585         float deltaphi = maxHBPhi - fPhi;
0586         if (fPhi > maxHBPhi) {
0587           deltaphi = fPhi - maxHBPhi;
0588         }
0589         if (deltaphi > PI) {
0590           deltaphi = 2.0 * PI - deltaphi;
0591         }
0592         float deltaeta = fEta - maxHBEta;
0593         Double_t r = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
0594
0595         HBCalREC.push_back(jhbhe->energy());
0596         HBCalR.push_back(r);
0597         HBCalSHE.push_back(fHBEnergySimHits[cell.rawId()]);
0598       }
0599
0600       if (cell.subdet() == sdHcalEC) {
0601         ++iHE;
0602
0603         const HcalGeometry* cellGeometry =
0604             dynamic_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, cell.subdet()));
0605         double fEta = cellGeometry->getPosition(cell).eta();
0606         double fPhi = cellGeometry->getPosition(cell).phi();
0607
0608         float deltaphi = maxHEPhi - fPhi;
0609         if (fPhi > maxHEPhi) {
0610           deltaphi = fPhi - maxHEPhi;
0611         }
0612         if (deltaphi > PI) {
0613           deltaphi = 2.0 * PI - deltaphi;
0614         }
0615         float deltaeta = fEta - maxHEEta;
0616         Double_t r = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
0617
0618         HECalREC.push_back(jhbhe->energy());
0619         HECalR.push_back(r);
0620         HECalSHE.push_back(fHEEnergySimHits[cell.rawId()]);
0621       }
0622     }
0623   }  // end loop through collection
0624
0625   if (verbosity > 1) {
0626     eventout += "\n          Number of HBRecHits collected:............ ";
0627     eventout += iHB;
0628   }
0629
0630   if (verbosity > 1) {
0631     eventout += "\n          Number of HERecHits collected:............ ";
0632     eventout += iHE;
0633   }
0634
0635   ////////////////////////
0636   // get HF information
0637   ///////////////////////
0638   std::vector<edm::Handle<HFRecHitCollection>> hf;
0639   iEvent.getManyByType(hf);
0640   if (!hf[0].isValid()) {
0641     edm::LogWarning(MsgLoggerCat) << "Unable to find any HFRecHitCollections in event!";
0642     return;
0643   }
0644   std::vector<edm::Handle<HFRecHitCollection>>::iterator ihf;
0645
0646   int iHF = 0;
0647   for (ihf = hf.begin(); ihf != hf.end(); ++ihf) {
0648     // find max values
0649     for (HFRecHitCollection::const_iterator jhf = (*ihf)->begin(); jhf != (*ihf)->end(); ++jhf) {
0650       HcalDetId cell(jhf->id());
0651
0652       if (cell.subdet() == sdHcalFwd) {
0653         auto cellGeometry = geometry->getSubdetectorGeometry(cell)->getGeometry(cell);
0654         double fEta = cellGeometry->getPosition().eta();
0655         double fPhi = cellGeometry->getPosition().phi();
0656         if ((jhf->energy()) > maxHFEnergy) {
0657           maxHFEnergy = jhf->energy();
0658           maxHFPhi = fPhi;
0659           maxHFEta = fEta;
0660         }
0661       }
0662     }  // end find max values
0663
0664     for (HFRecHitCollection::const_iterator jhf = (*ihf)->begin(); jhf != (*ihf)->end(); ++jhf) {
0665       HcalDetId cell(jhf->id());
0666
0667       if (cell.subdet() == sdHcalFwd) {
0668         ++iHF;
0669
0670         auto cellGeometry = geometry->getSubdetectorGeometry(cell)->getGeometry(cell);
0671         double fEta = cellGeometry->getPosition().eta();
0672         double fPhi = cellGeometry->getPosition().phi();
0673
0674         float deltaphi = maxHBPhi - fPhi;
0675         if (fPhi > maxHFPhi) {
0676           deltaphi = fPhi - maxHFPhi;
0677         }
0678         if (deltaphi > PI) {
0679           deltaphi = 2.0 * PI - deltaphi;
0680         }
0681         float deltaeta = fEta - maxHFEta;
0682         Double_t r = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
0683
0684         HFCalREC.push_back(jhf->energy());
0685         HFCalR.push_back(r);
0686         HFCalSHE.push_back(fHFEnergySimHits[cell.rawId()]);
0687       }
0688     }
0689   }  // end loop through collection
0690
0691   if (verbosity > 1) {
0692     eventout += "\n          Number of HFDigis collected:.............. ";
0693     eventout += iHF;
0694   }
0695
0696   ////////////////////////
0697   // get HO information
0698   ///////////////////////
0699   std::vector<edm::Handle<HORecHitCollection>> ho;
0700   iEvent.getManyByType(ho);
0701   if (!ho[0].isValid()) {
0702     edm::LogWarning(MsgLoggerCat) << "Unable to find any HORecHitCollections in event!";
0703     return;
0704   }
0705   std::vector<edm::Handle<HORecHitCollection>>::iterator iho;
0706
0707   int iHO = 0;
0708   for (iho = ho.begin(); iho != ho.end(); ++iho) {
0709     for (HORecHitCollection::const_iterator jho = (*iho)->begin(); jho != (*iho)->end(); ++jho) {
0710       HcalDetId cell(jho->id());
0711
0712       if (cell.subdet() == sdHcalOut) {
0713         ++iHO;
0714
0715         auto cellGeometry = geometry->getSubdetectorGeometry(cell)->getGeometry(cell);
0716         double fEta = cellGeometry->getPosition().eta();
0717         double fPhi = cellGeometry->getPosition().phi();
0718
0719         float deltaphi = maxHOPhi - fPhi;
0720         if (fPhi > maxHOPhi) {
0721           deltaphi = fPhi - maxHOPhi;
0722         }
0723         if (deltaphi > PI) {
0724           deltaphi = 2.0 * PI - deltaphi;
0725         }
0726         float deltaeta = fEta - maxHOEta;
0727         Double_t r = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
0728
0729         HOCalREC.push_back(jho->energy());
0730         HOCalR.push_back(r);
0731         HOCalSHE.push_back(fHOEnergySimHits[cell.rawId()]);
0732       }
0733     }
0734   }  // end loop through collection
0735
0736   if (verbosity > 1) {
0737     eventout += "\n          Number of HODigis collected:.............. ";
0738     eventout += iHO;
0739   }
0740
0741   if (verbosity > 0)
0742     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0743
0744   return;
0745 }
0746
0747 void GlobalRecHitsProducer::storeHCal(PGlobalRecHit& product) {
0748   std::string MsgLoggerCat = "GlobalRecHitsProducer_storeHCal";
0749
0750   if (verbosity > 2) {
0751     TString eventout("\n         nHBRecHits     = ");
0752     eventout += HBCalREC.size();
0753     for (unsigned int i = 0; i < HBCalREC.size(); ++i) {
0754       eventout += "\n      (REC, R, SHE) = (";
0755       eventout += HBCalREC[i];
0756       eventout += ", ";
0757       eventout += HBCalR[i];
0758       eventout += ", ";
0759       eventout += HBCalSHE[i];
0760       eventout += ")";
0761     }
0762     eventout += "\n         nHERecHits     = ";
0763     eventout += HECalREC.size();
0764     for (unsigned int i = 0; i < HECalREC.size(); ++i) {
0765       eventout += "\n      (REC, R, SHE) = (";
0766       eventout += HECalREC[i];
0767       eventout += ", ";
0768       eventout += HECalR[i];
0769       eventout += ", ";
0770       eventout += HECalSHE[i];
0771       eventout += ")";
0772     }
0773     eventout += "\n         nHFRecHits     = ";
0774     eventout += HFCalREC.size();
0775     for (unsigned int i = 0; i < HFCalREC.size(); ++i) {
0776       eventout += "\n      (REC, R, SHE) = (";
0777       eventout += HFCalREC[i];
0778       eventout += ", ";
0779       eventout += HFCalR[i];
0780       eventout += ", ";
0781       eventout += HFCalSHE[i];
0782       eventout += ")";
0783     }
0784     eventout += "\n         nHORecHits     = ";
0785     eventout += HOCalREC.size();
0786     for (unsigned int i = 0; i < HOCalREC.size(); ++i) {
0787       eventout += "\n      (REC, R, SHE) = (";
0788       eventout += HOCalREC[i];
0789       eventout += ", ";
0790       eventout += HOCalR[i];
0791       eventout += ", ";
0792       eventout += HOCalSHE[i];
0793       eventout += ")";
0794     }
0795
0796     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0797   }
0798
0799   product.putHBCalRecHits(HBCalREC, HBCalR, HBCalSHE);
0800   product.putHECalRecHits(HECalREC, HECalR, HECalSHE);
0801   product.putHOCalRecHits(HOCalREC, HOCalR, HOCalSHE);
0802   product.putHFCalRecHits(HFCalREC, HFCalR, HFCalSHE);
0803
0804   return;
0805 }
0806
0807 void GlobalRecHitsProducer::fillTrk(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0808   //Retrieve tracker topology from geometry
0809   const TrackerTopology* const tTopo = &iSetup.getData(tTopoToken_);
0810   ;
0811   std::string MsgLoggerCat = "GlobalRecHitsProducer_fillTrk";
0812   TString eventout;
0813   if (verbosity > 0)
0814     eventout = "\nGathering info:";
0815
0816   // get strip information
0817   edm::Handle<SiStripMatchedRecHit2DCollection> rechitsmatched;
0818   iEvent.getByToken(SiStripSrc_Token_, rechitsmatched);
0819   if (!rechitsmatched.isValid()) {
0820     edm::LogWarning(MsgLoggerCat) << "Unable to find stripmatchedrechits in event!";
0821     return;
0822   }
0823
0824   TrackerHitAssociator associate(iEvent, trackerHitAssociatorConfig_);
0825
0826   const auto& tGeomHandle = iSetup.getHandle(tGeomToken_);
0827   if (!tGeomHandle.isValid()) {
0828     edm::LogWarning(MsgLoggerCat) << "Unable to find TrackerDigiGeometry in event!";
0829     return;
0830   }
0831   const TrackerGeometry& tracker(*tGeomHandle);
0832
0833   int nStripBrl = 0, nStripFwd = 0;
0834
0835   // loop over det units
0836   for (TrackerGeometry::DetContainer::const_iterator it = tGeomHandle->dets().begin(); it != tGeomHandle->dets().end();
0837        ++it) {
0838     uint32_t myid = ((*it)->geographicalId()).rawId();
0839     DetId detid = ((*it)->geographicalId());
0840
0841     //loop over rechits-matched in the same subdetector
0842     SiStripMatchedRecHit2DCollection::const_iterator rechitmatchedMatch = rechitsmatched->find(detid);
0843
0844     if (rechitmatchedMatch != rechitsmatched->end()) {
0845       SiStripMatchedRecHit2DCollection::DetSet rechitmatchedRange = *rechitmatchedMatch;
0846       SiStripMatchedRecHit2DCollection::DetSet::const_iterator rechitmatchedRangeIteratorBegin =
0847           rechitmatchedRange.begin();
0848       SiStripMatchedRecHit2DCollection::DetSet::const_iterator rechitmatchedRangeIteratorEnd = rechitmatchedRange.end();
0849       SiStripMatchedRecHit2DCollection::DetSet::const_iterator itermatched;
0850
0851       for (itermatched = rechitmatchedRangeIteratorBegin; itermatched != rechitmatchedRangeIteratorEnd; ++itermatched) {
0852         SiStripMatchedRecHit2D const rechit = *itermatched;
0853         LocalPoint position = rechit.localPosition();
0854
0855         float mindist = 999999.;
0856         float distx = 999999.;
0857         float disty = 999999.;
0858         float dist = 999999.;
0859         std::pair<LocalPoint, LocalVector> closestPair;
0860         matched.clear();
0861
0862         float rechitmatchedx = position.x();
0863         float rechitmatchedy = position.y();
0864
0865         matched = associate.associateHit(rechit);
0866
0867         if (!matched.empty()) {
0868           //project simhit;
0869           const GluedGeomDet* gluedDet = (const GluedGeomDet*)tracker.idToDet(rechit.geographicalId());
0870           const StripGeomDetUnit* partnerstripdet = (StripGeomDetUnit*)gluedDet->stereoDet();
0871           std::pair<LocalPoint, LocalVector> hitPair;
0872
0873           for (std::vector<PSimHit>::const_iterator m = matched.begin(); m != matched.end(); m++) {
0874             //project simhit;
0875             hitPair = projectHit((*m), partnerstripdet, gluedDet->surface());
0876             distx = fabs(rechitmatchedx - hitPair.first.x());
0877             disty = fabs(rechitmatchedy - hitPair.first.y());
0878             dist = sqrt(distx * distx + disty * disty);
0879
0880             if (dist < mindist) {
0881               mindist = dist;
0882               closestPair = hitPair;
0883             }
0884           }
0885
0886           // get TIB
0887           if (detid.subdetId() == sdSiTIB) {
0888             ++nStripBrl;
0889
0890             if (tTopo->tibLayer(myid) == 1) {
0891               TIBL1RX.push_back(rechitmatchedx);
0892               TIBL1RY.push_back(rechitmatchedy);
0893               TIBL1SX.push_back(closestPair.first.x());
0894               TIBL1SY.push_back(closestPair.first.y());
0895             }
0896             if (tTopo->tibLayer(myid) == 2) {
0897               TIBL2RX.push_back(rechitmatchedx);
0898               TIBL2RY.push_back(rechitmatchedy);
0899               TIBL2SX.push_back(closestPair.first.x());
0900               TIBL2SY.push_back(closestPair.first.y());
0901             }
0902             if (tTopo->tibLayer(myid) == 3) {
0903               TIBL3RX.push_back(rechitmatchedx);
0904               TIBL3RY.push_back(rechitmatchedy);
0905               TIBL3SX.push_back(closestPair.first.x());
0906               TIBL3SY.push_back(closestPair.first.y());
0907             }
0908             if (tTopo->tibLayer(myid) == 4) {
0909               TIBL4RX.push_back(rechitmatchedx);
0910               TIBL4RY.push_back(rechitmatchedy);
0911               TIBL4SX.push_back(closestPair.first.x());
0912               TIBL4SY.push_back(closestPair.first.y());
0913             }
0914           }
0915
0916           // get TOB
0917           if (detid.subdetId() == sdSiTOB) {
0918             ++nStripBrl;
0919
0920             if (tTopo->tobLayer(myid) == 1) {
0921               TOBL1RX.push_back(rechitmatchedx);
0922               TOBL1RY.push_back(rechitmatchedy);
0923               TOBL1SX.push_back(closestPair.first.x());
0924               TOBL1SY.push_back(closestPair.first.y());
0925             }
0926             if (tTopo->tobLayer(myid) == 2) {
0927               TOBL2RX.push_back(rechitmatchedx);
0928               TOBL2RY.push_back(rechitmatchedy);
0929               TOBL2SX.push_back(closestPair.first.x());
0930               TOBL2SY.push_back(closestPair.first.y());
0931             }
0932             if (tTopo->tobLayer(myid) == 3) {
0933               TOBL3RX.push_back(rechitmatchedx);
0934               TOBL3RY.push_back(rechitmatchedy);
0935               TOBL3SX.push_back(closestPair.first.x());
0936               TOBL3SY.push_back(closestPair.first.y());
0937             }
0938             if (tTopo->tobLayer(myid) == 4) {
0939               TOBL4RX.push_back(rechitmatchedx);
0940               TOBL4RY.push_back(rechitmatchedy);
0941               TOBL4SX.push_back(closestPair.first.x());
0942               TOBL4SY.push_back(closestPair.first.y());
0943             }
0944           }
0945
0946           // get TID
0947           if (detid.subdetId() == sdSiTID) {
0948             ++nStripFwd;
0949
0950             if (tTopo->tidWheel(myid) == 1) {
0951               TIDW1RX.push_back(rechitmatchedx);
0952               TIDW1RY.push_back(rechitmatchedy);
0953               TIDW1SX.push_back(closestPair.first.x());
0954               TIDW1SY.push_back(closestPair.first.y());
0955             }
0956             if (tTopo->tidWheel(myid) == 2) {
0957               TIDW2RX.push_back(rechitmatchedx);
0958               TIDW2RY.push_back(rechitmatchedy);
0959               TIDW2SX.push_back(closestPair.first.x());
0960               TIDW2SY.push_back(closestPair.first.y());
0961             }
0962             if (tTopo->tidWheel(myid) == 3) {
0963               TIDW3RX.push_back(rechitmatchedx);
0964               TIDW3RY.push_back(rechitmatchedy);
0965               TIDW3SX.push_back(closestPair.first.x());
0966               TIDW3SY.push_back(closestPair.first.y());
0967             }
0968           }
0969
0970           // get TEC
0971           if (detid.subdetId() == sdSiTEC) {
0972             ++nStripFwd;
0973
0974             if (tTopo->tecWheel(myid) == 1) {
0975               TECW1RX.push_back(rechitmatchedx);
0976               TECW1RY.push_back(rechitmatchedy);
0977               TECW1SX.push_back(closestPair.first.x());
0978               TECW1SY.push_back(closestPair.first.y());
0979             }
0980             if (tTopo->tecWheel(myid) == 2) {
0981               TECW2RX.push_back(rechitmatchedx);
0982               TECW2RY.push_back(rechitmatchedy);
0983               TECW2SX.push_back(closestPair.first.x());
0984               TECW2SY.push_back(closestPair.first.y());
0985             }
0986             if (tTopo->tecWheel(myid) == 3) {
0987               TECW3RX.push_back(rechitmatchedx);
0988               TECW3RY.push_back(rechitmatchedy);
0989               TECW3SX.push_back(closestPair.first.x());
0990               TECW3SY.push_back(closestPair.first.y());
0991             }
0992             if (tTopo->tecWheel(myid) == 4) {
0993               TECW4RX.push_back(rechitmatchedx);
0994               TECW4RY.push_back(rechitmatchedy);
0995               TECW4SX.push_back(closestPair.first.x());
0996               TECW4SY.push_back(closestPair.first.y());
0997             }
0998             if (tTopo->tecWheel(myid) == 5) {
0999               TECW5RX.push_back(rechitmatchedx);
1000               TECW5RY.push_back(rechitmatchedy);
1001               TECW5SX.push_back(closestPair.first.x());
1002               TECW5SY.push_back(closestPair.first.y());
1003             }
1004             if (tTopo->tecWheel(myid) == 6) {
1005               TECW6RX.push_back(rechitmatchedx);
1006               TECW6RY.push_back(rechitmatchedy);
1007               TECW6SX.push_back(closestPair.first.x());
1008               TECW6SY.push_back(closestPair.first.y());
1009             }
1010             if (tTopo->tecWheel(myid) == 7) {
1011               TECW7RX.push_back(rechitmatchedx);
1012               TECW7RY.push_back(rechitmatchedy);
1013               TECW7SX.push_back(closestPair.first.x());
1014               TECW7SY.push_back(closestPair.first.y());
1015             }
1016             if (tTopo->tecWheel(myid) == 8) {
1017               TECW8RX.push_back(rechitmatchedx);
1018               TECW8RY.push_back(rechitmatchedy);
1019               TECW8SX.push_back(closestPair.first.x());
1020               TECW8SY.push_back(closestPair.first.y());
1021             }
1022           }
1023
1024         }  // end if matched empty
1025       }
1026     }
1027   }  // end loop over det units
1028
1029   if (verbosity > 1) {
1030     eventout += "\n          Number of BrlStripRecHits collected:...... ";
1031     eventout += nStripBrl;
1032   }
1033
1034   if (verbosity > 1) {
1035     eventout += "\n          Number of FrwdStripRecHits collected:..... ";
1036     eventout += nStripFwd;
1037   }
1038
1039   // get pixel information
1040   //Get RecHits
1041   edm::Handle<SiPixelRecHitCollection> recHitColl;
1042   iEvent.getByToken(SiPxlSrc_Token_, recHitColl);
1043   if (!recHitColl.isValid()) {
1044     edm::LogWarning(MsgLoggerCat) << "Unable to find SiPixelRecHitCollection in event!";
1045     return;
1046   }
1047
1048   int nPxlBrl = 0, nPxlFwd = 0;
1049   //iterate over detunits
1050   for (TrackerGeometry::DetContainer::const_iterator it = tGeomHandle->dets().begin(); it != tGeomHandle->dets().end();
1051        ++it) {
1052     uint32_t myid = ((*it)->geographicalId()).rawId();
1053     DetId detId = ((*it)->geographicalId());
1054     int subid = detId.subdetId();
1055
1056     if (!((subid == sdPxlBrl) || (subid == sdPxlFwd)))
1057       continue;
1058
1059     //const PixelGeomDetUnit * theGeomDet =
1060     //  dynamic_cast<const PixelGeomDetUnit*>(theTracker.idToDet(detId) );
1061
1062     SiPixelRecHitCollection::const_iterator pixeldet = recHitColl->find(detId);
1063     if (pixeldet == recHitColl->end())
1064       continue;
1065     SiPixelRecHitCollection::DetSet pixelrechitRange = *pixeldet;
1066     SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorBegin = pixelrechitRange.begin();
1067     SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorEnd = pixelrechitRange.end();
1068     SiPixelRecHitCollection::DetSet::const_iterator pixeliter = pixelrechitRangeIteratorBegin;
1069     std::vector<PSimHit> matched;
1070
1071     //----Loop over rechits for this detId
1072     for (; pixeliter != pixelrechitRangeIteratorEnd; ++pixeliter) {
1073       matched.clear();
1074       matched = associate.associateHit(*pixeliter);
1075
1076       if (!matched.empty()) {
1077         float closest = 9999.9;
1078         //std::vector<PSimHit>::const_iterator closestit = matched.begin();
1079         LocalPoint lp = pixeliter->localPosition();
1080         float rechit_x = lp.x();
1081         float rechit_y = lp.y();
1082
1083         float sim_x = 0.;
1084         float sim_y = 0.;
1085
1086         //loop over sim hits and fill closet
1087         for (std::vector<PSimHit>::const_iterator m = matched.begin(); m != matched.end(); ++m) {
1088           float sim_x1 = (*m).entryPoint().x();
1089           float sim_x2 = (*m).exitPoint().x();
1090           float sim_xpos = 0.5 * (sim_x1 + sim_x2);
1091
1092           float sim_y1 = (*m).entryPoint().y();
1093           float sim_y2 = (*m).exitPoint().y();
1094           float sim_ypos = 0.5 * (sim_y1 + sim_y2);
1095
1096           float x_res = fabs(sim_xpos - rechit_x);
1097           float y_res = fabs(sim_ypos - rechit_y);
1098
1099           float dist = sqrt(x_res * x_res + y_res * y_res);
1100
1101           if (dist < closest) {
1102             closest = dist;
1103             sim_x = sim_xpos;
1104             sim_y = sim_ypos;
1105           }
1106         }  // end sim hit loop
1107
1108         // get Barrel pixels
1109         if (subid == sdPxlBrl) {
1110           ++nPxlBrl;
1111
1112           if (tTopo->pxbLayer(myid) == 1) {
1113             BRL1RX.push_back(rechit_x);
1114             BRL1RY.push_back(rechit_y);
1115             BRL1SX.push_back(sim_x);
1116             BRL1SY.push_back(sim_y);
1117           }
1118           if (tTopo->pxbLayer(myid) == 2) {
1119             BRL2RX.push_back(rechit_x);
1120             BRL2RY.push_back(rechit_y);
1121             BRL2SX.push_back(sim_x);
1122             BRL2SY.push_back(sim_y);
1123           }
1124           if (tTopo->pxbLayer(myid) == 3) {
1125             BRL3RX.push_back(rechit_x);
1126             BRL3RY.push_back(rechit_y);
1127             BRL3SX.push_back(sim_x);
1128             BRL3SY.push_back(sim_y);
1129           }
1130         }
1131
1132         // get Forward pixels
1133         if (subid == sdPxlFwd) {
1134           ++nPxlFwd;
1135
1136           if (tTopo->pxfDisk(myid) == 1) {
1137             if (tTopo->pxfSide(myid) == 1) {
1138               FWD1nRX.push_back(rechit_x);
1139               FWD1nRY.push_back(rechit_y);
1140               FWD1nSX.push_back(sim_x);
1141               FWD1nSY.push_back(sim_y);
1142             }
1143             if (tTopo->pxfSide(myid) == 2) {
1144               FWD1pRX.push_back(rechit_x);
1145               FWD1pRY.push_back(rechit_y);
1146               FWD1pSX.push_back(sim_x);
1147               FWD1pSY.push_back(sim_y);
1148             }
1149           }
1150           if (tTopo->pxfDisk(myid) == 2) {
1151             if (tTopo->pxfSide(myid) == 1) {
1152               FWD2nRX.push_back(rechit_x);
1153               FWD2nRY.push_back(rechit_y);
1154               FWD2nSX.push_back(sim_x);
1155               FWD2nSY.push_back(sim_y);
1156             }
1157             if (tTopo->pxfSide(myid) == 2) {
1158               FWD2pRX.push_back(rechit_x);
1159               FWD2pRY.push_back(rechit_y);
1160               FWD2pSX.push_back(sim_x);
1161               FWD2pSY.push_back(sim_y);
1162             }
1163           }
1164         }
1165       }  // end matched emtpy
1166     }    // <-----end rechit loop
1167   }      // <------ end detunit loop
1168
1169   if (verbosity > 1) {
1170     eventout += "\n          Number of BrlPixelRecHits collected:...... ";
1171     eventout += nPxlBrl;
1172   }
1173
1174   if (verbosity > 1) {
1175     eventout += "\n          Number of FrwdPixelRecHits collected:..... ";
1176     eventout += nPxlFwd;
1177   }
1178
1179   if (verbosity > 0)
1180     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1181
1182   return;
1183 }
1184
1185 void GlobalRecHitsProducer::storeTrk(PGlobalRecHit& product) {
1186   std::string MsgLoggerCat = "GlobalRecHitsProducer_storeTrk";
1187
1188   if (verbosity > 2) {
1189     // strip output
1190     TString eventout("\n         nTIBL1     = ");
1191     eventout += TIBL1RX.size();
1192     for (unsigned int i = 0; i < TIBL1RX.size(); ++i) {
1193       eventout += "\n      (RX, RY, SX, SY) = (";
1194       eventout += TIBL1RX[i];
1195       eventout += ", ";
1196       eventout += TIBL1RY[i];
1197       eventout += ", ";
1198       eventout += TIBL1SX[i];
1199       eventout += ", ";
1200       eventout += TIBL1SY[i];
1201       eventout += ")";
1202     }
1203     eventout += "\n         nTIBL2     = ";
1204     eventout += TIBL2RX.size();
1205     for (unsigned int i = 0; i < TIBL2RX.size(); ++i) {
1206       eventout += "\n      (RX, RY, SX, SY) = (";
1207       eventout += TIBL2RX[i];
1208       eventout += ", ";
1209       eventout += TIBL2RY[i];
1210       eventout += ", ";
1211       eventout += TIBL2SX[i];
1212       eventout += ", ";
1213       eventout += TIBL2SY[i];
1214       eventout += ")";
1215     }
1216     eventout += "\n         nTIBL3     = ";
1217     eventout += TIBL3RX.size();
1218     for (unsigned int i = 0; i < TIBL3RX.size(); ++i) {
1219       eventout += "\n      (RX, RY, SX, SY) = (";
1220       eventout += TIBL3RX[i];
1221       eventout += ", ";
1222       eventout += TIBL3RY[i];
1223       eventout += ", ";
1224       eventout += TIBL3SX[i];
1225       eventout += ", ";
1226       eventout += TIBL3SY[i];
1227       eventout += ")";
1228     }
1229     eventout += "\n         nTIBL4     = ";
1230     eventout += TIBL4RX.size();
1231     for (unsigned int i = 0; i < TIBL4RX.size(); ++i) {
1232       eventout += "\n      (RX, RY, SX, SY) = (";
1233       eventout += TIBL4RX[i];
1234       eventout += ", ";
1235       eventout += TIBL4RY[i];
1236       eventout += ", ";
1237       eventout += TIBL4SX[i];
1238       eventout += ", ";
1239       eventout += TIBL4SY[i];
1240       eventout += ")";
1241     }
1242     eventout += "\n         nTOBL1     = ";
1243     eventout += TOBL1RX.size();
1244     for (unsigned int i = 0; i < TOBL1RX.size(); ++i) {
1245       eventout += "\n      (RX, RY, SX, SY) = (";
1246       eventout += TOBL1RX[i];
1247       eventout += ", ";
1248       eventout += TOBL1RY[i];
1249       eventout += ", ";
1250       eventout += TOBL1SX[i];
1251       eventout += ", ";
1252       eventout += TOBL1SY[i];
1253       eventout += ")";
1254     }
1255     eventout += "\n         nTOBL2     = ";
1256     eventout += TOBL2RX.size();
1257     for (unsigned int i = 0; i < TOBL2RX.size(); ++i) {
1258       eventout += "\n      (RX, RY, SX, SY) = (";
1259       eventout += TOBL2RX[i];
1260       eventout += ", ";
1261       eventout += TOBL2RY[i];
1262       eventout += ", ";
1263       eventout += TOBL2SX[i];
1264       eventout += ", ";
1265       eventout += TOBL2SY[i];
1266       eventout += ")";
1267     }
1268     eventout += "\n         nTOBL3     = ";
1269     eventout += TOBL3RX.size();
1270     for (unsigned int i = 0; i < TOBL3RX.size(); ++i) {
1271       eventout += "\n      (RX, RY, SX, SY) = (";
1272       eventout += TOBL3RX[i];
1273       eventout += ", ";
1274       eventout += TOBL3RY[i];
1275       eventout += ", ";
1276       eventout += TOBL3SX[i];
1277       eventout += ", ";
1278       eventout += TOBL3SY[i];
1279       eventout += ")";
1280     }
1281     eventout += "\n         nTOBL4     = ";
1282     eventout += TOBL4RX.size();
1283     for (unsigned int i = 0; i < TOBL4RX.size(); ++i) {
1284       eventout += "\n      (RX, RY, SX, SY) = (";
1285       eventout += TOBL4RX[i];
1286       eventout += ", ";
1287       eventout += TOBL4RY[i];
1288       eventout += ", ";
1289       eventout += TOBL4SX[i];
1290       eventout += ", ";
1291       eventout += TOBL4SY[i];
1292       eventout += ")";
1293     }
1294     eventout += "\n         nTIDW1     = ";
1295     eventout += TIDW1RX.size();
1296     for (unsigned int i = 0; i < TIDW1RX.size(); ++i) {
1297       eventout += "\n      (RX, RY, SX, SY) = (";
1298       eventout += TIDW1RX[i];
1299       eventout += ", ";
1300       eventout += TIDW1RY[i];
1301       eventout += ", ";
1302       eventout += TIDW1SX[i];
1303       eventout += ", ";
1304       eventout += TIDW1SY[i];
1305       eventout += ")";
1306     }
1307     eventout += "\n         nTIDW2     = ";
1308     eventout += TIDW2RX.size();
1309     for (unsigned int i = 0; i < TIDW2RX.size(); ++i) {
1310       eventout += "\n      (RX, RY, SX, SY) = (";
1311       eventout += TIDW2RX[i];
1312       eventout += ", ";
1313       eventout += TIDW2RY[i];
1314       eventout += ", ";
1315       eventout += TIDW2SX[i];
1316       eventout += ", ";
1317       eventout += TIDW2SY[i];
1318       eventout += ")";
1319     }
1320     eventout += "\n         nTIDW3     = ";
1321     eventout += TIDW3RX.size();
1322     for (unsigned int i = 0; i < TIDW3RX.size(); ++i) {
1323       eventout += "\n      (RX, RY, SX, SY) = (";
1324       eventout += TIDW3RX[i];
1325       eventout += ", ";
1326       eventout += TIDW3RY[i];
1327       eventout += ", ";
1328       eventout += TIDW3SX[i];
1329       eventout += ", ";
1330       eventout += TIDW3SY[i];
1331       eventout += ")";
1332     }
1333     eventout += "\n         nTECW1     = ";
1334     eventout += TECW1RX.size();
1335     for (unsigned int i = 0; i < TECW1RX.size(); ++i) {
1336       eventout += "\n      (RX, RY, SX, SY) = (";
1337       eventout += TECW1RX[i];
1338       eventout += ", ";
1339       eventout += TECW1RY[i];
1340       eventout += ", ";
1341       eventout += TECW1SX[i];
1342       eventout += ", ";
1343       eventout += TECW1SY[i];
1344       eventout += ")";
1345     }
1346     eventout += "\n         nTECW2     = ";
1347     eventout += TECW2RX.size();
1348     for (unsigned int i = 0; i < TECW2RX.size(); ++i) {
1349       eventout += "\n      (RX, RY, SX, SY) = (";
1350       eventout += TECW2RX[i];
1351       eventout += ", ";
1352       eventout += TECW2RY[i];
1353       eventout += ", ";
1354       eventout += TECW2SX[i];
1355       eventout += ", ";
1356       eventout += TECW2SY[i];
1357       eventout += ")";
1358     }
1359     eventout += "\n         nTECW3     = ";
1360     eventout += TECW3RX.size();
1361     for (unsigned int i = 0; i < TECW3RX.size(); ++i) {
1362       eventout += "\n      (RX, RY, SX, SY) = (";
1363       eventout += TECW3RX[i];
1364       eventout += ", ";
1365       eventout += TECW3RY[i];
1366       eventout += ", ";
1367       eventout += TECW3SX[i];
1368       eventout += ", ";
1369       eventout += TECW3SY[i];
1370       eventout += ")";
1371     }
1372     eventout += "\n         nTECW4     = ";
1373     eventout += TECW4RX.size();
1374     for (unsigned int i = 0; i < TECW4RX.size(); ++i) {
1375       eventout += "\n      (RX, RY, SX, SY) = (";
1376       eventout += TECW4RX[i];
1377       eventout += ", ";
1378       eventout += TECW4RY[i];
1379       eventout += ", ";
1380       eventout += TECW4SX[i];
1381       eventout += ", ";
1382       eventout += TECW4SY[i];
1383       eventout += ")";
1384     }
1385     eventout += "\n         nTECW5     = ";
1386     eventout += TECW5RX.size();
1387     for (unsigned int i = 0; i < TECW5RX.size(); ++i) {
1388       eventout += "\n      (RX, RY, SX, SY) = (";
1389       eventout += TECW5RX[i];
1390       eventout += ", ";
1391       eventout += TECW5RY[i];
1392       eventout += ", ";
1393       eventout += TECW5SX[i];
1394       eventout += ", ";
1395       eventout += TECW5SY[i];
1396       eventout += ")";
1397     }
1398     eventout += "\n         nTECW6     = ";
1399     eventout += TECW6RX.size();
1400     for (unsigned int i = 0; i < TECW6RX.size(); ++i) {
1401       eventout += "\n      (RX, RY, SX, SY) = (";
1402       eventout += TECW6RX[i];
1403       eventout += ", ";
1404       eventout += TECW6RY[i];
1405       eventout += ", ";
1406       eventout += TECW6SX[i];
1407       eventout += ", ";
1408       eventout += TECW6SY[i];
1409       eventout += ")";
1410     }
1411     eventout += "\n         nTECW7     = ";
1412     eventout += TECW7RX.size();
1413     for (unsigned int i = 0; i < TECW7RX.size(); ++i) {
1414       eventout += "\n      (RX, RY, SX, SY) = (";
1415       eventout += TECW7RX[i];
1416       eventout += ", ";
1417       eventout += TECW7RY[i];
1418       eventout += ", ";
1419       eventout += TECW7SX[i];
1420       eventout += ", ";
1421       eventout += TECW7SY[i];
1422       eventout += ")";
1423     }
1424     eventout += "\n         nTECW8     = ";
1425     eventout += TECW8RX.size();
1426     for (unsigned int i = 0; i < TECW8RX.size(); ++i) {
1427       eventout += "\n      (RX, RY, SX, SY) = (";
1428       eventout += TECW8RX[i];
1429       eventout += ", ";
1430       eventout += TECW8RY[i];
1431       eventout += ", ";
1432       eventout += TECW8SX[i];
1433       eventout += ", ";
1434       eventout += TECW8SY[i];
1435       eventout += ")";
1436     }
1437
1438     // pixel output
1439     eventout += "\n         nBRL1     = ";
1440     eventout += BRL1RX.size();
1441     for (unsigned int i = 0; i < BRL1RX.size(); ++i) {
1442       eventout += "\n      (RX, RY, SX, SY) = (";
1443       eventout += BRL1RX[i];
1444       eventout += ", ";
1445       eventout += BRL1RY[i];
1446       eventout += ", ";
1447       eventout += BRL1SX[i];
1448       eventout += ", ";
1449       eventout += BRL1SY[i];
1450       eventout += ")";
1451     }
1452     eventout += "\n         nBRL2     = ";
1453     eventout += BRL2RX.size();
1454     for (unsigned int i = 0; i < BRL2RX.size(); ++i) {
1455       eventout += "\n      (RX, RY, SX, SY) = (";
1456       eventout += BRL2RX[i];
1457       eventout += ", ";
1458       eventout += BRL2RY[i];
1459       eventout += ", ";
1460       eventout += BRL2SX[i];
1461       eventout += ", ";
1462       eventout += BRL2SY[i];
1463       eventout += ")";
1464     }
1465     eventout += "\n         nBRL3     = ";
1466     eventout += BRL3RX.size();
1467     for (unsigned int i = 0; i < BRL3RX.size(); ++i) {
1468       eventout += "\n      (RX, RY, SX, SY) = (";
1469       eventout += BRL3RX[i];
1470       eventout += ", ";
1471       eventout += BRL3RY[i];
1472       eventout += ", ";
1473       eventout += BRL3SX[i];
1474       eventout += ", ";
1475       eventout += BRL3SY[i];
1476       eventout += ")";
1477     }
1478     eventout += "\n         nFWD1p     = ";
1479     eventout += FWD1pRX.size();
1480     for (unsigned int i = 0; i < FWD1pRX.size(); ++i) {
1481       eventout += "\n      (RX, RY, SX, SY) = (";
1482       eventout += FWD1pRX[i];
1483       eventout += ", ";
1484       eventout += FWD1pRY[i];
1485       eventout += ", ";
1486       eventout += FWD1pSX[i];
1487       eventout += ", ";
1488       eventout += FWD1pSY[i];
1489       eventout += ")";
1490     }
1491     eventout += "\n         nFWD1n     = ";
1492     eventout += FWD1nRX.size();
1493     for (unsigned int i = 0; i < FWD1nRX.size(); ++i) {
1494       eventout += "\n      (RX, RY, SX, SY) = (";
1495       eventout += FWD1nRX[i];
1496       eventout += ", ";
1497       eventout += FWD1nRY[i];
1498       eventout += ", ";
1499       eventout += FWD1nSX[i];
1500       eventout += ", ";
1501       eventout += FWD1nSY[i];
1502       eventout += ")";
1503     }
1504     eventout += "\n         nFWD2p     = ";
1505     eventout += FWD2pRX.size();
1506     for (unsigned int i = 0; i < FWD2pRX.size(); ++i) {
1507       eventout += "\n      (RX, RY, SX, SY) = (";
1508       eventout += FWD2pRX[i];
1509       eventout += ", ";
1510       eventout += FWD2pRY[i];
1511       eventout += ", ";
1512       eventout += FWD2pSX[i];
1513       eventout += ", ";
1514       eventout += FWD2pSY[i];
1515       eventout += ")";
1516     }
1517     eventout += "\n         nFWD2p     = ";
1518     eventout += FWD2nRX.size();
1519     for (unsigned int i = 0; i < FWD2nRX.size(); ++i) {
1520       eventout += "\n      (RX, RY, SX, SY) = (";
1521       eventout += FWD2nRX[i];
1522       eventout += ", ";
1523       eventout += FWD2nRY[i];
1524       eventout += ", ";
1525       eventout += FWD2nSX[i];
1526       eventout += ", ";
1527       eventout += FWD2nSY[i];
1528       eventout += ")";
1529     }
1530
1531     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1532   }
1533
1534   // strip output
1535   product.putTIBL1RecHits(TIBL1RX, TIBL1RY, TIBL1SX, TIBL1SY);
1536   product.putTIBL2RecHits(TIBL2RX, TIBL2RY, TIBL2SX, TIBL2SY);
1537   product.putTIBL3RecHits(TIBL3RX, TIBL3RY, TIBL3SX, TIBL3SY);
1538   product.putTIBL4RecHits(TIBL4RX, TIBL4RY, TIBL4SX, TIBL4SY);
1539   product.putTOBL1RecHits(TOBL1RX, TOBL1RY, TOBL1SX, TOBL1SY);
1540   product.putTOBL2RecHits(TOBL2RX, TOBL2RY, TOBL2SX, TOBL2SY);
1541   product.putTOBL3RecHits(TOBL3RX, TOBL3RY, TOBL3SX, TOBL3SY);
1542   product.putTOBL4RecHits(TOBL4RX, TOBL4RY, TOBL4SX, TOBL4SY);
1543   product.putTIDW1RecHits(TIDW1RX, TIDW1RY, TIDW1SX, TIDW1SY);
1544   product.putTIDW2RecHits(TIDW2RX, TIDW2RY, TIDW2SX, TIDW2SY);
1545   product.putTIDW3RecHits(TIDW3RX, TIDW3RY, TIDW3SX, TIDW3SY);
1546   product.putTECW1RecHits(TECW1RX, TECW1RY, TECW1SX, TECW1SY);
1547   product.putTECW2RecHits(TECW2RX, TECW2RY, TECW2SX, TECW2SY);
1548   product.putTECW3RecHits(TECW3RX, TECW3RY, TECW3SX, TECW3SY);
1549   product.putTECW4RecHits(TECW4RX, TECW4RY, TECW4SX, TECW4SY);
1550   product.putTECW5RecHits(TECW5RX, TECW5RY, TECW5SX, TECW5SY);
1551   product.putTECW6RecHits(TECW6RX, TECW6RY, TECW6SX, TECW6SY);
1552   product.putTECW7RecHits(TECW7RX, TECW7RY, TECW7SX, TECW7SY);
1553   product.putTECW8RecHits(TECW8RX, TECW8RY, TECW8SX, TECW8SY);
1554
1555   // pixel output
1556   product.putBRL1RecHits(BRL1RX, BRL1RY, BRL1SX, BRL1SY);
1557   product.putBRL2RecHits(BRL2RX, BRL2RY, BRL2SX, BRL2SY);
1558   product.putBRL3RecHits(BRL3RX, BRL3RY, BRL3SX, BRL3SY);
1559   product.putFWD1pRecHits(FWD1pRX, FWD1pRY, FWD1pSX, FWD1pSY);
1560   product.putFWD1nRecHits(FWD1nRX, FWD1nRY, FWD1nSX, FWD1nSY);
1561   product.putFWD2pRecHits(FWD2pRX, FWD2pRY, FWD2pSX, FWD2pSY);
1562   product.putFWD2nRecHits(FWD2nRX, FWD2nRY, FWD2nSX, FWD2nSY);
1563
1564   return;
1565 }
1566
1567 void GlobalRecHitsProducer::fillMuon(edm::Event& iEvent, const edm::EventSetup& iSetup) {
1568   std::string MsgLoggerCat = "GlobalRecHitsProducer_fillMuon";
1569
1570   TString eventout;
1571   if (verbosity > 0)
1572     eventout = "\nGathering info:";
1573
1574   // get DT information
1575   const auto& dtGeom = iSetup.getHandle(dtGeomToken_);
1576   if (!dtGeom.isValid()) {
1577     edm::LogWarning(MsgLoggerCat) << "Unable to find DTMuonGeometryRecord in event!";
1578     return;
1579   }
1580
1581   edm::Handle<edm::PSimHitContainer> dtsimHits;
1582   iEvent.getByToken(MuDTSimSrc_Token_, dtsimHits);
1583   if (!dtsimHits.isValid()) {
1584     edm::LogWarning(MsgLoggerCat) << "Unable to find dtsimHits in event!";
1585     return;
1586   }
1587
1588   std::map<DTWireId, edm::PSimHitContainer> simHitsPerWire =
1589       DTHitQualityUtils::mapSimHitsPerWire(*(dtsimHits.product()));
1590
1591   edm::Handle<DTRecHitCollection> dtRecHits;
1592   iEvent.getByToken(MuDTSrc_Token_, dtRecHits);
1593   if (!dtRecHits.isValid()) {
1594     edm::LogWarning(MsgLoggerCat) << "Unable to find dtRecHits in event!";
1595     return;
1596   }
1597
1598   std::map<DTWireId, std::vector<DTRecHit1DPair>> recHitsPerWire = map1DRecHitsPerWire(dtRecHits.product());
1599
1600   int nDt = compute(dtGeom.product(), simHitsPerWire, recHitsPerWire, 1);
1601
1602   if (verbosity > 1) {
1603     eventout += "\n          Number of DtMuonRecHits collected:........ ";
1604     eventout += nDt;
1605   }
1606
1607   // get CSC Strip information
1608   // get map of sim hits
1609   theMap.clear();
1610   //edm::Handle<CrossingFrame> cf;
1611   edm::Handle<CrossingFrame<PSimHit>> cf;
1612   //iEvent.getByType(cf);
1613   //if (!cf.isValid()) {
1614   //  edm::LogWarning(MsgLoggerCat)
1615   //    << "Unable to find CrossingFrame in event!";
1616   //  return;
1617   //}
1618   //MixCollection<PSimHit> simHits(cf.product(), "MuonCSCHits");
1619   iEvent.getByToken(MuCSCHits_Token_, cf);
1620   if (!cf.isValid()) {
1621     edm::LogWarning(MsgLoggerCat) << "Unable to find muo CSC  crossingFrame in event!";
1622     return;
1623   }
1624   MixCollection<PSimHit> simHits(cf.product());
1625
1626   // arrange the hits by detUnit
1627   for (MixCollection<PSimHit>::MixItr hitItr = simHits.begin(); hitItr != simHits.end(); ++hitItr) {
1628     theMap[hitItr->detUnitId()].push_back(*hitItr);
1629   }
1630
1631   // get geometry
1632   const auto& hGeom = iSetup.getHandle(cscGeomToken_);
1633   if (!hGeom.isValid()) {
1634     edm::LogWarning(MsgLoggerCat) << "Unable to find CSCMuonGeometryRecord in event!";
1635     return;
1636   }
1637   const CSCGeometry* theCSCGeometry = &*hGeom;
1638
1639   // get rechits
1640   edm::Handle<CSCRecHit2DCollection> hRecHits;
1641   iEvent.getByToken(MuCSCSrc_Token_, hRecHits);
1642   if (!hRecHits.isValid()) {
1643     edm::LogWarning(MsgLoggerCat) << "Unable to find CSC RecHits in event!";
1644     return;
1645   }
1646   const CSCRecHit2DCollection* cscRecHits = hRecHits.product();
1647
1648   int nCSC = 0;
1649   for (CSCRecHit2DCollection::const_iterator recHitItr = cscRecHits->begin(); recHitItr != cscRecHits->end();
1650        ++recHitItr) {
1651     int detId = (*recHitItr).cscDetId().rawId();
1652
1653     edm::PSimHitContainer simHits;
1654     std::map<int, edm::PSimHitContainer>::const_iterator mapItr = theMap.find(detId);
1655     if (mapItr != theMap.end()) {
1656       simHits = mapItr->second;
1657     }
1658
1659     if (simHits.size() == 1) {
1660       ++nCSC;
1661
1662       const GeomDetUnit* detUnit = theCSCGeometry->idToDetUnit(CSCDetId(detId));
1663       const CSCLayer* layer = dynamic_cast<const CSCLayer*>(detUnit);
1664
1665       int chamberType = layer->chamber()->specs()->chamberType();
1666       plotResolution(simHits[0], *recHitItr, layer, chamberType);
1667     }
1668   }
1669
1670   if (verbosity > 1) {
1671     eventout += "\n          Number of CSCRecHits collected:........... ";
1672     eventout += nCSC;
1673   }
1674
1675   // get RPC information
1676   std::map<double, int> mapsim, maprec;
1677   std::map<int, double> nmapsim, nmaprec;
1678
1679   const auto& rpcGeom = iSetup.getHandle(rpcGeomToken_);
1680   if (!rpcGeom.isValid()) {
1681     edm::LogWarning(MsgLoggerCat) << "Unable to find RPCMuonGeometryRecord in event!";
1682     return;
1683   }
1684
1685   edm::Handle<edm::PSimHitContainer> simHit;
1686   iEvent.getByToken(MuRPCSimSrc_Token_, simHit);
1687   if (!simHit.isValid()) {
1688     edm::LogWarning(MsgLoggerCat) << "Unable to find RPCSimHit in event!";
1689     return;
1690   }
1691
1692   edm::Handle<RPCRecHitCollection> recHit;
1693   iEvent.getByToken(MuRPCSrc_Token_, recHit);
1694   if (!simHit.isValid()) {
1695     edm::LogWarning(MsgLoggerCat) << "Unable to find RPCRecHit in event!";
1696     return;
1697   }
1698
1699   int nRPC = 0;
1700   RPCRecHitCollection::const_iterator recIt;
1701   int nrec = 0;
1702   for (recIt = recHit->begin(); recIt != recHit->end(); ++recIt) {
1703     RPCDetId Rid = (RPCDetId)(*recIt).rpcId();
1704     const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(Rid));
1705     if (roll->isForward()) {
1706       if (verbosity > 1) {
1707         eventout += "\n          Number of RPCRecHits collected:........... ";
1708         eventout += nRPC;
1709       }
1710
1711       if (verbosity > 0)
1712         edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1713       return;
1714     }
1715     nrec = nrec + 1;
1716     LocalPoint rhitlocal = (*recIt).localPosition();
1717     double rhitlocalx = rhitlocal.x();
1718     maprec[rhitlocalx] = nrec;
1719   }
1720
1721   int i = 0;
1722   for (std::map<double, int>::iterator iter = maprec.begin(); iter != maprec.end(); ++iter) {
1723     i = i + 1;
1724     nmaprec[i] = (*iter).first;
1725   }
1726
1727   edm::PSimHitContainer::const_iterator simIt;
1728   int nsim = 0;
1729   for (simIt = simHit->begin(); simIt != simHit->end(); simIt++) {
1730     int ptype = (*simIt).particleType();
1731     //RPCDetId Rsid = (RPCDetId)(*simIt).detUnitId();
1732     if (ptype == 13 || ptype == -13) {
1733       nsim = nsim + 1;
1734       LocalPoint shitlocal = (*simIt).localPosition();
1735       double shitlocalx = shitlocal.x();
1736       mapsim[shitlocalx] = nsim;
1737     }
1738   }
1739
1740   i = 0;
1741   for (std::map<double, int>::iterator iter = mapsim.begin(); iter != mapsim.end(); ++iter) {
1742     i = i + 1;
1743     nmapsim[i] = (*iter).first;
1744   }
1745
1746   if (nsim == nrec) {
1747     for (int r = 0; r < nsim; r++) {
1748       ++nRPC;
1749       RPCRHX.push_back(nmaprec[r + 1]);
1750       RPCSHX.push_back(nmapsim[r + 1]);
1751     }
1752   }
1753
1754   if (verbosity > 1) {
1755     eventout += "\n          Number of RPCRecHits collected:........... ";
1756     eventout += nRPC;
1757   }
1758
1759   if (verbosity > 0)
1760     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1761
1762   return;
1763 }
1764
1765 void GlobalRecHitsProducer::storeMuon(PGlobalRecHit& product) {
1766   std::string MsgLoggerCat = "GlobalRecHitsProducer_storeMuon";
1767
1768   if (verbosity > 2) {
1769     // dt output
1770     TString eventout("\n         nDT     = ");
1771     eventout += DTRHD.size();
1772     for (unsigned int i = 0; i < DTRHD.size(); ++i) {
1773       eventout += "\n      (RHD, SHD) = (";
1774       eventout += DTRHD[i];
1775       eventout += ", ";
1776       eventout += DTSHD[i];
1777       eventout += ")";
1778     }
1779
1780     // CSC Strip
1781     eventout += "\n         nCSC     = ";
1782     eventout += CSCRHPHI.size();
1783     for (unsigned int i = 0; i < CSCRHPHI.size(); ++i) {
1784       eventout += "\n      (rhphi, rhperp, shphi) = (";
1785       eventout += CSCRHPHI[i];
1786       eventout += ", ";
1787       eventout += CSCRHPERP[i];
1788       eventout += ", ";
1789       eventout += CSCSHPHI[i];
1790       eventout += ")";
1791     }
1792
1793     // RPC
1794     eventout += "\n         nRPC     = ";
1795     eventout += RPCRHX.size();
1796     for (unsigned int i = 0; i < RPCRHX.size(); ++i) {
1797       eventout += "\n      (rhx, shx) = (";
1798       eventout += RPCRHX[i];
1799       eventout += ", ";
1800       eventout += RPCSHX[i];
1801       eventout += ")";
1802     }
1803
1804     edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1805   }
1806
1807   product.putDTRecHits(DTRHD, DTSHD);
1808
1809   product.putCSCRecHits(CSCRHPHI, CSCRHPERP, CSCSHPHI);
1810
1811   product.putRPCRecHits(RPCRHX, RPCSHX);
1812
1813   return;
1814 }
1815
1816 void GlobalRecHitsProducer::clear() {
1817   std::string MsgLoggerCat = "GlobalRecHitsProducer_clear";
1818
1819   if (verbosity > 0)
1820     edm::LogInfo(MsgLoggerCat) << "Clearing event holders";
1821
1822   // reset electromagnetic info
1823   // EE info
1824   EERE.clear();
1825   EESHE.clear();
1826   // EB info
1827   EBRE.clear();
1828   EBSHE.clear();
1829   // ES info
1830   ESRE.clear();
1831   ESSHE.clear();
1832
1833   // reset HCal Info
1834   HBCalREC.clear();
1835   HBCalR.clear();
1836   HBCalSHE.clear();
1837   HECalREC.clear();
1838   HECalR.clear();
1839   HECalSHE.clear();
1840   HOCalREC.clear();
1841   HOCalR.clear();
1842   HOCalSHE.clear();
1843   HFCalREC.clear();
1844   HFCalR.clear();
1845   HFCalSHE.clear();
1846
1847   // reset Track Info
1848   TIBL1RX.clear();
1849   TIBL2RX.clear();
1850   TIBL3RX.clear();
1851   TIBL4RX.clear();
1852   TIBL1RY.clear();
1853   TIBL2RY.clear();
1854   TIBL3RY.clear();
1855   TIBL4RY.clear();
1856   TIBL1SX.clear();
1857   TIBL2SX.clear();
1858   TIBL3SX.clear();
1859   TIBL4SX.clear();
1860   TIBL1SY.clear();
1861   TIBL2SY.clear();
1862   TIBL3SY.clear();
1863   TIBL4SY.clear();
1864
1865   TOBL1RX.clear();
1866   TOBL2RX.clear();
1867   TOBL3RX.clear();
1868   TOBL4RX.clear();
1869   TOBL1RY.clear();
1870   TOBL2RY.clear();
1871   TOBL3RY.clear();
1872   TOBL4RY.clear();
1873   TOBL1SX.clear();
1874   TOBL2SX.clear();
1875   TOBL3SX.clear();
1876   TOBL4SX.clear();
1877   TOBL1SY.clear();
1878   TOBL2SY.clear();
1879   TOBL3SY.clear();
1880   TOBL4SY.clear();
1881
1882   TIDW1RX.clear();
1883   TIDW2RX.clear();
1884   TIDW3RX.clear();
1885   TIDW1RY.clear();
1886   TIDW2RY.clear();
1887   TIDW3RY.clear();
1888   TIDW1SX.clear();
1889   TIDW2SX.clear();
1890   TIDW3SX.clear();
1891   TIDW1SY.clear();
1892   TIDW2SY.clear();
1893   TIDW3SY.clear();
1894
1895   TECW1RX.clear();
1896   TECW2RX.clear();
1897   TECW3RX.clear();
1898   TECW4RX.clear();
1899   TECW5RX.clear();
1900   TECW6RX.clear();
1901   TECW7RX.clear();
1902   TECW8RX.clear();
1903   TECW1RY.clear();
1904   TECW2RY.clear();
1905   TECW3RY.clear();
1906   TECW4RY.clear();
1907   TECW5RY.clear();
1908   TECW6RY.clear();
1909   TECW7RY.clear();
1910   TECW8RY.clear();
1911   TECW1SX.clear();
1912   TECW2SX.clear();
1913   TECW3SX.clear();
1914   TECW4SX.clear();
1915   TECW5SX.clear();
1916   TECW6SX.clear();
1917   TECW7SX.clear();
1918   TECW8SX.clear();
1919   TECW1SY.clear();
1920   TECW2SY.clear();
1921   TECW3SY.clear();
1922   TECW4SY.clear();
1923   TECW5SY.clear();
1924   TECW6SY.clear();
1925   TECW7SY.clear();
1926   TECW8SY.clear();
1927
1928   BRL1RX.clear();
1929   BRL1RY.clear();
1930   BRL1SX.clear();
1931   BRL1SY.clear();
1932   BRL2RX.clear();
1933   BRL2RY.clear();
1934   BRL2SX.clear();
1935   BRL2SY.clear();
1936   BRL3RX.clear();
1937   BRL3RY.clear();
1938   BRL3SX.clear();
1939   BRL3SY.clear();
1940
1941   FWD1pRX.clear();
1942   FWD1pRY.clear();
1943   FWD1pSX.clear();
1944   FWD1pSY.clear();
1945   FWD1nRX.clear();
1946   FWD1nRY.clear();
1947   FWD1nSX.clear();
1948   FWD1nSY.clear();
1949   FWD2pRX.clear();
1950   FWD2pRY.clear();
1951   FWD2pSX.clear();
1952   FWD2pSY.clear();
1953   FWD2nRX.clear();
1954   FWD2nRY.clear();
1955   FWD2nSX.clear();
1956   FWD2nSY.clear();
1957
1958   //muon clear
1959   DTRHD.clear();
1960   DTSHD.clear();
1961
1962   CSCRHPHI.clear();
1963   CSCRHPERP.clear();
1964   CSCSHPHI.clear();
1965
1966   RPCRHX.clear();
1967   RPCSHX.clear();
1968
1969   return;
1970 }
1971
1972 //needed by to do the residual for matched hits in SiStrip
1973 std::pair<LocalPoint, LocalVector> GlobalRecHitsProducer::projectHit(const PSimHit& hit,
1974                                                                      const StripGeomDetUnit* stripDet,
1975                                                                      const BoundPlane& plane) {
1976   const StripTopology& topol = stripDet->specificTopology();
1977   GlobalPoint globalpos = stripDet->surface().toGlobal(hit.localPosition());
1978   LocalPoint localHit = plane.toLocal(globalpos);
1979   //track direction
1980   LocalVector locdir = hit.localDirection();
1981   //rotate track in new frame
1982
1983   GlobalVector globaldir = stripDet->surface().toGlobal(locdir);
1984   LocalVector dir = plane.toLocal(globaldir);
1985   float scale = -localHit.z() / dir.z();
1986
1987   LocalPoint projectedPos = localHit + scale * dir;
1988
1989   float selfAngle = topol.stripAngle(topol.strip(hit.localPosition()));
1990
1991   // vector along strip in hit frame
1992   LocalVector stripDir(sin(selfAngle), cos(selfAngle), 0);
1993
1994   LocalVector localStripDir(plane.toLocal(stripDet->surface().toGlobal(stripDir)));
1995
1996   return std::pair<LocalPoint, LocalVector>(projectedPos, localStripDir);
1997 }
1998
1999 // Return a map between DTRecHit1DPair and wireId
2000 std::map<DTWireId, std::vector<DTRecHit1DPair>> GlobalRecHitsProducer::map1DRecHitsPerWire(
2001     const DTRecHitCollection* dt1DRecHitPairs) {
2002   std::map<DTWireId, std::vector<DTRecHit1DPair>> ret;
2003
2004   for (DTRecHitCollection::const_iterator rechit = dt1DRecHitPairs->begin(); rechit != dt1DRecHitPairs->end();
2005        rechit++) {
2006     ret[(*rechit).wireId()].push_back(*rechit);
2007   }
2008
2009   return ret;
2010 }
2011
2012 // Compute SimHit distance from wire (cm)
2013 float GlobalRecHitsProducer::simHitDistFromWire(const DTLayer* layer, DTWireId wireId, const PSimHit& hit) {
2014   float xwire = layer->specificTopology().wirePosition(wireId.wire());
2015   LocalPoint entryP = hit.entryPoint();
2016   LocalPoint exitP = hit.exitPoint();
2017   float xEntry = entryP.x() - xwire;
2018   float xExit = exitP.x() - xwire;
2019
2020   //FIXME: check...
2021   return fabs(xEntry - (entryP.z() * (xExit - xEntry)) / (exitP.z() - entryP.z()));
2022 }
2023
2024 // Find the RecHit closest to the muon SimHit
2025 template <typename type>
2026 const type* GlobalRecHitsProducer::findBestRecHit(const DTLayer* layer,
2027                                                   DTWireId wireId,
2028                                                   const std::vector<type>& recHits,
2029                                                   const float simHitDist) {
2030   float res = 99999;
2031   const type* theBestRecHit = nullptr;
2032   // Loop over RecHits within the cell
2033   for (typename std::vector<type>::const_iterator recHit = recHits.begin(); recHit != recHits.end(); recHit++) {
2034     float distTmp = recHitDistFromWire(*recHit, layer);
2035     if (fabs(distTmp - simHitDist) < res) {
2036       res = fabs(distTmp - simHitDist);
2037       theBestRecHit = &(*recHit);
2038     }
2039   }  // End of loop over RecHits within the cell
2040
2041   return theBestRecHit;
2042 }
2043
2044 // Compute the distance from wire (cm) of a hits in a DTRecHit1DPair
2045 float GlobalRecHitsProducer::recHitDistFromWire(const DTRecHit1DPair& hitPair, const DTLayer* layer) {
2046   // Compute the rechit distance from wire
2047   return fabs(hitPair.localPosition(DTEnums::Left).x() - hitPair.localPosition(DTEnums::Right).x()) / 2.;
2048 }
2049
2050 // Compute the distance from wire (cm) of a hits in a DTRecHit1D
2051 float GlobalRecHitsProducer::recHitDistFromWire(const DTRecHit1D& recHit, const DTLayer* layer) {
2052   return fabs(recHit.localPosition().x() - layer->specificTopology().wirePosition(recHit.wireId().wire()));
2053 }
2054
2055 template <typename type>
2056 int GlobalRecHitsProducer::compute(const DTGeometry* dtGeom,
2057                                    const std::map<DTWireId, std::vector<PSimHit>>& _simHitsPerWire,
2058                                    const std::map<DTWireId, std::vector<type>>& _recHitsPerWire,
2059                                    int step) {
2060   std::map<DTWireId, std::vector<PSimHit>> simHitsPerWire = _simHitsPerWire;
2061   std::map<DTWireId, std::vector<type>> recHitsPerWire = _recHitsPerWire;
2062   int nDt = 0;
2063   // Loop over cells with a muon SimHit
2064   for (std::map<DTWireId, std::vector<PSimHit>>::const_iterator wireAndSHits = simHitsPerWire.begin();
2065        wireAndSHits != simHitsPerWire.end();
2066        wireAndSHits++) {
2067     DTWireId wireId = (*wireAndSHits).first;
2068     std::vector<PSimHit> simHitsInCell = (*wireAndSHits).second;
2069
2070     // Get the layer
2071     const DTLayer* layer = dtGeom->layer(wireId);
2072
2073     // Look for a mu hit in the cell
2074     const PSimHit* muSimHit = DTHitQualityUtils::findMuSimHit(simHitsInCell);
2075     if (muSimHit == nullptr) {
2076       continue;  // Skip this cell
2077     }
2078
2079     // Find the distance of the simhit from the wire
2080     float simHitWireDist = simHitDistFromWire(layer, wireId, *muSimHit);
2081     // Skip simhits out of the cell
2082     if (simHitWireDist > 2.1) {
2083       continue;  // Skip this cell
2084     }
2085     //GlobalPoint simHitGlobalPos = layer->toGlobal(muSimHit->localPosition());
2086
2087     // Look for RecHits in the same cell
2088     if (recHitsPerWire.find(wireId) == recHitsPerWire.end()) {
2089       continue;  // No RecHit found in this cell
2090     } else {
2091       // vector<type> recHits = (*wireAndRecHits).second;
2092       std::vector<type> recHits = recHitsPerWire[wireId];
2093
2094       // Find the best RecHit
2095       const type* theBestRecHit = findBestRecHit(layer, wireId, recHits, simHitWireDist);
2096
2097       float recHitWireDist = recHitDistFromWire(*theBestRecHit, layer);
2098
2099       ++nDt;
2100
2101       DTRHD.push_back(recHitWireDist);
2102       DTSHD.push_back(simHitWireDist);
2103
2104     }  // find rechits
2105   }    // loop over simhits
2106
2107   return nDt;
2108 }
2109
2110 void GlobalRecHitsProducer::plotResolution(const PSimHit& simHit,
2111                                            const CSCRecHit2D& recHit,
2112                                            const CSCLayer* layer,
2113                                            int chamberType) {
2114   GlobalPoint simHitPos = layer->toGlobal(simHit.localPosition());
2115   GlobalPoint recHitPos = layer->toGlobal(recHit.localPosition());
2116
2117   CSCRHPHI.push_back(recHitPos.phi());
2118   CSCRHPERP.push_back(recHitPos.perp());
2119   CSCSHPHI.push_back(simHitPos.phi());
2120 }
2121
2122 //define this as a plug-in
2123 //DEFINE_FWK_MODULE(GlobalRecHitsProducer);