File indexing completed on 2023-03-17 11:27:39
0001
0002
0003
0004
0005
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
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
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
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
0082
0083 verbosity %= 10;
0084
0085
0086 produces<PGlobalRecHit>(label);
0087
0088
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
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
0142 ++count;
0143
0144
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
0158 clear();
0159
0160
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
0175 eventout += AllProv[i]->moduleLabel();
0176 eventout += "\n ProductID : ";
0177
0178 eventout += AllProv[i]->productID().id();
0179 eventout += "\n ClassName : ";
0180
0181 eventout += AllProv[i]->className();
0182 eventout += "\n InstanceName : ";
0183
0184 eventout += AllProv[i]->productInstanceName();
0185 eventout += "\n BranchName : ";
0186
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
0197
0198 fillECal(iEvent, iSetup);
0199
0200 fillHCal(iEvent, iSetup);
0201
0202 fillTrk(iEvent, iSetup);
0203
0204 fillMuon(iEvent, iSetup);
0205
0206 if (verbosity > 0)
0207 edm::LogInfo(MsgLoggerCat) << "Done gathering data from event.";
0208
0209
0210 std::unique_ptr<PGlobalRecHit> pOut(new PGlobalRecHit);
0211
0212 if (verbosity > 2)
0213 edm::LogInfo(MsgLoggerCat) << "Saving event contents:";
0214
0215
0216
0217 storeECal(*pOut);
0218
0219 storeHCal(*pOut);
0220
0221 storeTrk(*pOut);
0222
0223 storeMuon(*pOut);
0224
0225
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
0239
0240 edm::Handle<CrossingFrame<PCaloHit>> crossingFrame;
0241
0242
0243
0244
0245
0246
0247
0248
0249
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
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
0272
0273
0274 std::unique_ptr<MixCollection<PCaloHit>> barrelHits(new MixCollection<PCaloHit>(crossingFrame.product()));
0275
0276
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
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
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
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
0333
0334
0335 std::unique_ptr<MixCollection<PCaloHit>> endcapHits(new MixCollection<PCaloHit>(crossingFrame.product()));
0336
0337
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
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
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
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
0387
0388
0389 std::unique_ptr<MixCollection<PCaloHit>> preshowerHits(new MixCollection<PCaloHit>(crossingFrame.product()));
0390
0391
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
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
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
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
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
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
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 }
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 }
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
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
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 }
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 }
0690
0691 if (verbosity > 1) {
0692 eventout += "\n Number of HFDigis collected:.............. ";
0693 eventout += iHF;
0694 }
0695
0696
0697
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 }
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
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
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
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
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
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
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
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
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
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
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 }
1025 }
1026 }
1027 }
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
1040
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
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
1060
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
1072 for (; pixeliter != pixelrechitRangeIteratorEnd; ++pixeliter) {
1073 matched.clear();
1074 matched = associate.associateHit(*pixeliter);
1075
1076 if (!matched.empty()) {
1077 float closest = 9999.9;
1078
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
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 }
1107
1108
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
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 }
1166 }
1167 }
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
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
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
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
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
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
1608
1609 theMap.clear();
1610
1611 edm::Handle<CrossingFrame<PSimHit>> cf;
1612
1613
1614
1615
1616
1617
1618
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
1627 for (MixCollection<PSimHit>::MixItr hitItr = simHits.begin(); hitItr != simHits.end(); ++hitItr) {
1628 theMap[hitItr->detUnitId()].push_back(*hitItr);
1629 }
1630
1631
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
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
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
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
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
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
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
1823
1824 EERE.clear();
1825 EESHE.clear();
1826
1827 EBRE.clear();
1828 EBSHE.clear();
1829
1830 ESRE.clear();
1831 ESSHE.clear();
1832
1833
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
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
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
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
1980 LocalVector locdir = hit.localDirection();
1981
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
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
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
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
2021 return fabs(xEntry - (entryP.z() * (xExit - xEntry)) / (exitP.z() - entryP.z()));
2022 }
2023
2024
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
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 }
2040
2041 return theBestRecHit;
2042 }
2043
2044
2045 float GlobalRecHitsProducer::recHitDistFromWire(const DTRecHit1DPair& hitPair, const DTLayer* layer) {
2046
2047 return fabs(hitPair.localPosition(DTEnums::Left).x() - hitPair.localPosition(DTEnums::Right).x()) / 2.;
2048 }
2049
2050
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
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
2071 const DTLayer* layer = dtGeom->layer(wireId);
2072
2073
2074 const PSimHit* muSimHit = DTHitQualityUtils::findMuSimHit(simHitsInCell);
2075 if (muSimHit == nullptr) {
2076 continue;
2077 }
2078
2079
2080 float simHitWireDist = simHitDistFromWire(layer, wireId, *muSimHit);
2081
2082 if (simHitWireDist > 2.1) {
2083 continue;
2084 }
2085
2086
2087
2088 if (recHitsPerWire.find(wireId) == recHitsPerWire.end()) {
2089 continue;
2090 } else {
2091
2092 std::vector<type> recHits = recHitsPerWire[wireId];
2093
2094
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 }
2105 }
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
2123