Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:32:25 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!