File indexing completed on 2022-02-14 23:39:57
0001
0002
0003
0004
0005
0006
0007
0008 #include "Validation/GlobalRecHits/interface/GlobalRecHitsAnalyzer.h"
0009 #include "DQMServices/Core/interface/DQMStore.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 using namespace std;
0015
0016 GlobalRecHitsAnalyzer::GlobalRecHitsAnalyzer(const edm::ParameterSet& iPSet)
0017 : fName(""),
0018 verbosity(0),
0019 frequency(0),
0020 label(""),
0021 getAllProvenances(false),
0022 printProvenanceInfo(false),
0023 trackerHitAssociatorConfig_(iPSet, consumesCollector()),
0024 caloGeomToken_(esConsumes()),
0025 tTopoToken_(esConsumes()),
0026 tGeomToken_(esConsumes()),
0027 dtGeomToken_(esConsumes()),
0028 cscGeomToken_(esConsumes()),
0029 rpcGeomToken_(esConsumes()),
0030 count(0) {
0031 consumesMany<edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>>();
0032 consumesMany<edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>>();
0033 consumesMany<edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>>();
0034 std::string MsgLoggerCat = "GlobalRecHitsAnalyzer_GlobalRecHitsAnalyzer";
0035
0036
0037 fName = iPSet.getUntrackedParameter<std::string>("Name");
0038 verbosity = iPSet.getUntrackedParameter<int>("Verbosity");
0039 frequency = iPSet.getUntrackedParameter<int>("Frequency");
0040 edm::ParameterSet m_Prov = iPSet.getParameter<edm::ParameterSet>("ProvenanceLookup");
0041 getAllProvenances = m_Prov.getUntrackedParameter<bool>("GetAllProvenances");
0042 printProvenanceInfo = m_Prov.getUntrackedParameter<bool>("PrintProvenanceInfo");
0043 hitsProducer = iPSet.getParameter<std::string>("hitsProducer");
0044
0045
0046 ECalEBSrc_ = iPSet.getParameter<edm::InputTag>("ECalEBSrc");
0047 ECalUncalEBSrc_ = iPSet.getParameter<edm::InputTag>("ECalUncalEBSrc");
0048 ECalEESrc_ = iPSet.getParameter<edm::InputTag>("ECalEESrc");
0049 ECalUncalEESrc_ = iPSet.getParameter<edm::InputTag>("ECalUncalEESrc");
0050 ECalESSrc_ = iPSet.getParameter<edm::InputTag>("ECalESSrc");
0051 HCalSrc_ = iPSet.getParameter<edm::InputTag>("HCalSrc");
0052 SiStripSrc_ = iPSet.getParameter<edm::InputTag>("SiStripSrc");
0053 SiPxlSrc_ = iPSet.getParameter<edm::InputTag>("SiPxlSrc");
0054 MuDTSrc_ = iPSet.getParameter<edm::InputTag>("MuDTSrc");
0055 MuDTSimSrc_ = iPSet.getParameter<edm::InputTag>("MuDTSimSrc");
0056 MuCSCSrc_ = iPSet.getParameter<edm::InputTag>("MuCSCSrc");
0057 MuRPCSrc_ = iPSet.getParameter<edm::InputTag>("MuRPCSrc");
0058 MuRPCSimSrc_ = iPSet.getParameter<edm::InputTag>("MuRPCSimSrc");
0059
0060
0061 ECalUncalEBSrc_Token_ = consumes<EBUncalibratedRecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalUncalEBSrc"));
0062 ECalUncalEESrc_Token_ = consumes<EEUncalibratedRecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalUncalEESrc"));
0063 ECalEBSrc_Token_ = consumes<EBRecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalEBSrc"));
0064 ECalEESrc_Token_ = consumes<EERecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalEESrc"));
0065 ECalESSrc_Token_ = consumes<ESRecHitCollection>(iPSet.getParameter<edm::InputTag>("ECalESSrc"));
0066 HCalSrc_Token_ = consumes<edm::PCaloHitContainer>(iPSet.getParameter<edm::InputTag>("HCalSrc"));
0067 SiStripSrc_Token_ = consumes<SiStripMatchedRecHit2DCollection>(iPSet.getParameter<edm::InputTag>("SiStripSrc"));
0068 SiPxlSrc_Token_ = consumes<SiPixelRecHitCollection>(iPSet.getParameter<edm::InputTag>("SiPxlSrc"));
0069
0070 MuDTSrc_Token_ = consumes<DTRecHitCollection>(iPSet.getParameter<edm::InputTag>("MuDTSrc"));
0071 MuDTSimSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("MuDTSimSrc"));
0072
0073 MuCSCSrc_Token_ = consumes<CSCRecHit2DCollection>(iPSet.getParameter<edm::InputTag>("MuCSCSrc"));
0074 MuCSCHits_Token_ = consumes<CrossingFrame<PSimHit>>(
0075 edm::InputTag(std::string("mix"), iPSet.getParameter<std::string>("hitsProducer") + std::string("MuonCSCHits")));
0076
0077 MuRPCSrc_Token_ = consumes<RPCRecHitCollection>(iPSet.getParameter<edm::InputTag>("MuRPCSrc"));
0078 MuRPCSimSrc_Token_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("MuRPCSimSrc"));
0079
0080 EBHits_Token_ = consumes<CrossingFrame<PCaloHit>>(
0081 edm::InputTag(std::string("mix"), iPSet.getParameter<std::string>("hitsProducer") + std::string("EcalHitsEB")));
0082 EEHits_Token_ = consumes<CrossingFrame<PCaloHit>>(
0083 edm::InputTag(std::string("mix"), iPSet.getParameter<std::string>("hitsProducer") + std::string("EcalHitsEE")));
0084 ESHits_Token_ = consumes<CrossingFrame<PCaloHit>>(
0085 edm::InputTag(std::string("mix"), iPSet.getParameter<std::string>("hitsProducer") + std::string("EcalHitsES")));
0086
0087
0088
0089 verbosity %= 10;
0090
0091
0092
0093
0094
0095 if (verbosity >= 0) {
0096 edm::LogInfo(MsgLoggerCat) << "\n===============================\n"
0097 << "Initialized as EDProducer with parameter values:\n"
0098 << " Name = " << fName << "\n"
0099 << " Verbosity = " << verbosity << "\n"
0100 << " Frequency = " << frequency << "\n"
0101 << " GetProv = " << getAllProvenances << "\n"
0102 << " PrintProv = " << printProvenanceInfo << "\n"
0103 << " ECalEBSrc = " << ECalEBSrc_.label() << ":" << ECalEBSrc_.instance() << "\n"
0104 << " ECalUncalEBSrc = " << ECalUncalEBSrc_.label() << ":"
0105 << ECalUncalEBSrc_.instance() << "\n"
0106 << " ECalEESrc = " << ECalEESrc_.label() << ":" << ECalUncalEESrc_.instance()
0107 << "\n"
0108 << " ECalUncalEESrc = " << ECalUncalEESrc_.label() << ":" << ECalEESrc_.instance()
0109 << "\n"
0110 << " ECalESSrc = " << ECalESSrc_.label() << ":" << ECalESSrc_.instance() << "\n"
0111 << " HCalSrc = " << HCalSrc_.label() << ":" << HCalSrc_.instance() << "\n"
0112 << " SiStripSrc = " << SiStripSrc_.label() << ":" << SiStripSrc_.instance()
0113 << "\n"
0114 << " SiPixelSrc = " << SiPxlSrc_.label() << ":" << SiPxlSrc_.instance() << "\n"
0115 << " MuDTSrc = " << MuDTSrc_.label() << ":" << MuDTSrc_.instance() << "\n"
0116 << " MuDTSimSrc = " << MuDTSimSrc_.label() << ":" << MuDTSimSrc_.instance()
0117 << "\n"
0118 << " MuCSCSrc = " << MuCSCSrc_.label() << ":" << MuCSCSrc_.instance() << "\n"
0119 << " MuRPCSrc = " << MuRPCSrc_.label() << ":" << MuRPCSrc_.instance() << "\n"
0120 << " MuRPCSimSrc = " << MuRPCSimSrc_.label() << ":" << MuRPCSimSrc_.instance()
0121 << "\n"
0122 << "===============================\n";
0123 }
0124 }
0125
0126 GlobalRecHitsAnalyzer::~GlobalRecHitsAnalyzer() {}
0127
0128 void GlobalRecHitsAnalyzer::bookHistograms(DQMStore::IBooker& iBooker, edm::Run const&, edm::EventSetup const&) {
0129
0130 string SiStripString[19] = {"TECW1",
0131 "TECW2",
0132 "TECW3",
0133 "TECW4",
0134 "TECW5",
0135 "TECW6",
0136 "TECW7",
0137 "TECW8",
0138 "TIBL1",
0139 "TIBL2",
0140 "TIBL3",
0141 "TIBL4",
0142 "TIDW1",
0143 "TIDW2",
0144 "TIDW3",
0145 "TOBL1",
0146 "TOBL2",
0147 "TOBL3",
0148 "TOBL4"};
0149 for (int i = 0; i < 19; ++i) {
0150 mehSiStripn[i] = nullptr;
0151 mehSiStripResX[i] = nullptr;
0152 mehSiStripResY[i] = nullptr;
0153 }
0154 string hcharname, hchartitle;
0155 iBooker.setCurrentFolder("GlobalRecHitsV/SiStrips");
0156 for (int amend = 0; amend < 19; ++amend) {
0157 hcharname = "hSiStripn_" + SiStripString[amend];
0158 hchartitle = SiStripString[amend] + " rechits";
0159 mehSiStripn[amend] = iBooker.book1D(hcharname, hchartitle, 200, 0., 200.);
0160 mehSiStripn[amend]->setAxisTitle("Number of hits in " + SiStripString[amend], 1);
0161 mehSiStripn[amend]->setAxisTitle("Count", 2);
0162 hcharname = "hSiStripResX_" + SiStripString[amend];
0163 hchartitle = SiStripString[amend] + " rechit x resolution";
0164 mehSiStripResX[amend] = iBooker.book1D(hcharname, hchartitle, 200, -0.02, .02);
0165 mehSiStripResX[amend]->setAxisTitle("X-resolution in " + SiStripString[amend], 1);
0166 mehSiStripResX[amend]->setAxisTitle("Count", 2);
0167 hcharname = "hSiStripResY_" + SiStripString[amend];
0168 hchartitle = SiStripString[amend] + " rechit y resolution";
0169 mehSiStripResY[amend] = iBooker.book1D(hcharname, hchartitle, 200, -0.02, .02);
0170 mehSiStripResY[amend]->setAxisTitle("Y-resolution in " + SiStripString[amend], 1);
0171 mehSiStripResY[amend]->setAxisTitle("Count", 2);
0172 }
0173
0174
0175
0176 string HCalString[4] = {"HB", "HE", "HF", "HO"};
0177 float HCalnUpper[4] = {3000., 3000., 3000., 3000.};
0178 float HCalnLower[4] = {0., 0., 0., 0.};
0179 for (int j = 0; j < 4; ++j) {
0180 mehHcaln[j] = nullptr;
0181 mehHcalRes[j] = nullptr;
0182 }
0183
0184 iBooker.setCurrentFolder("GlobalRecHitsV/HCals");
0185 for (int amend = 0; amend < 4; ++amend) {
0186 hcharname = "hHcaln_" + HCalString[amend];
0187 hchartitle = HCalString[amend] + " rechits";
0188 mehHcaln[amend] = iBooker.book1D(hcharname, hchartitle, 1000, HCalnLower[amend], HCalnUpper[amend]);
0189 mehHcaln[amend]->setAxisTitle("Number of RecHits", 1);
0190 mehHcaln[amend]->setAxisTitle("Count", 2);
0191 hcharname = "hHcalRes_" + HCalString[amend];
0192 hchartitle = HCalString[amend] + " rechit resolution";
0193 mehHcalRes[amend] = iBooker.book1D(hcharname, hchartitle, 25, -2., 2.);
0194 mehHcalRes[amend]->setAxisTitle("RecHit E - SimHit E", 1);
0195 mehHcalRes[amend]->setAxisTitle("Count", 2);
0196 }
0197
0198
0199 string ECalString[3] = {"EB", "EE", "ES"};
0200 int ECalnBins[3] = {1000, 3000, 150};
0201 float ECalnUpper[3] = {20000., 62000., 3000.};
0202 float ECalnLower[3] = {0., 0., 0.};
0203 int ECalResBins[3] = {200, 200, 200};
0204 float ECalResUpper[3] = {1., 0.3, .0002};
0205 float ECalResLower[3] = {-1., -0.3, -.0002};
0206 for (int i = 0; i < 3; ++i) {
0207 mehEcaln[i] = nullptr;
0208 mehEcalRes[i] = nullptr;
0209 }
0210 iBooker.setCurrentFolder("GlobalRecHitsV/ECals");
0211
0212 for (int amend = 0; amend < 3; ++amend) {
0213 hcharname = "hEcaln_" + ECalString[amend];
0214 hchartitle = ECalString[amend] + " rechits";
0215 mehEcaln[amend] = iBooker.book1D(hcharname, hchartitle, ECalnBins[amend], ECalnLower[amend], ECalnUpper[amend]);
0216 mehEcaln[amend]->setAxisTitle("Number of RecHits", 1);
0217 mehEcaln[amend]->setAxisTitle("Count", 2);
0218 hcharname = "hEcalRes_" + ECalString[amend];
0219 hchartitle = ECalString[amend] + " rechit resolution";
0220 mehEcalRes[amend] =
0221 iBooker.book1D(hcharname, hchartitle, ECalResBins[amend], ECalResLower[amend], ECalResUpper[amend]);
0222 mehEcalRes[amend]->setAxisTitle("RecHit E - SimHit E", 1);
0223 mehEcalRes[amend]->setAxisTitle("Count", 2);
0224 }
0225
0226
0227 string SiPixelString[7] = {"BRL1", "BRL2", "BRL3", "FWD1n", "FWD1p", "FWD2n", "FWD2p"};
0228 for (int j = 0; j < 7; ++j) {
0229 mehSiPixeln[j] = nullptr;
0230 mehSiPixelResX[j] = nullptr;
0231 mehSiPixelResY[j] = nullptr;
0232 }
0233
0234 iBooker.setCurrentFolder("GlobalRecHitsV/SiPixels");
0235 for (int amend = 0; amend < 7; ++amend) {
0236 hcharname = "hSiPixeln_" + SiPixelString[amend];
0237 hchartitle = SiPixelString[amend] + " rechits";
0238 mehSiPixeln[amend] = iBooker.book1D(hcharname, hchartitle, 200, 0., 200.);
0239 mehSiPixeln[amend]->setAxisTitle("Number of hits in " + SiPixelString[amend], 1);
0240 mehSiPixeln[amend]->setAxisTitle("Count", 2);
0241 hcharname = "hSiPixelResX_" + SiPixelString[amend];
0242 hchartitle = SiPixelString[amend] + " rechit x resolution";
0243 mehSiPixelResX[amend] = iBooker.book1D(hcharname, hchartitle, 200, -0.02, .02);
0244 mehSiPixelResX[amend]->setAxisTitle("X-resolution in " + SiPixelString[amend], 1);
0245 mehSiPixelResX[amend]->setAxisTitle("Count", 2);
0246 hcharname = "hSiPixelResY_" + SiPixelString[amend];
0247 hchartitle = SiPixelString[amend] + " rechit y resolution";
0248
0249 mehSiPixelResY[amend] = iBooker.book1D(hcharname, hchartitle, 200, -0.02, .02);
0250 mehSiPixelResY[amend]->setAxisTitle("Y-resolution in " + SiPixelString[amend], 1);
0251 mehSiPixelResY[amend]->setAxisTitle("Count", 2);
0252 }
0253
0254
0255 iBooker.setCurrentFolder("GlobalRecHitsV/Muons");
0256
0257 mehDtMuonn = nullptr;
0258 mehCSCn = nullptr;
0259 mehRPCn = nullptr;
0260
0261 string n_List[3] = {"hDtMuonn", "hCSCn", "hRPCn"};
0262 string hist_string[3] = {"Dt", "CSC", "RPC"};
0263
0264 for (int amend = 0; amend < 3; ++amend) {
0265 hchartitle = hist_string[amend] + " rechits";
0266 if (amend == 0) {
0267 mehDtMuonn = iBooker.book1D(n_List[amend], hchartitle, 50, 0., 500.);
0268 mehDtMuonn->setAxisTitle("Number of Rechits", 1);
0269 mehDtMuonn->setAxisTitle("Count", 2);
0270 }
0271 if (amend == 1) {
0272 mehCSCn = iBooker.book1D(n_List[amend], hchartitle, 50, 0., 500.);
0273 mehCSCn->setAxisTitle("Number of Rechits", 1);
0274 mehCSCn->setAxisTitle("Count", 2);
0275 }
0276 if (amend == 2) {
0277 mehRPCn = iBooker.book1D(n_List[amend], hchartitle, 50, 0., 500.);
0278 mehRPCn->setAxisTitle("Number of Rechits", 1);
0279 mehRPCn->setAxisTitle("Count", 2);
0280 }
0281 }
0282
0283 mehDtMuonRes = nullptr;
0284 mehCSCResRDPhi = nullptr;
0285 mehRPCResX = nullptr;
0286
0287 hcharname = "hDtMuonRes";
0288 hchartitle = "DT wire distance resolution";
0289 mehDtMuonRes = iBooker.book1D(hcharname, hchartitle, 200, -0.2, 0.2);
0290 hcharname = "CSCResRDPhi";
0291 hchartitle = "CSC perp*dphi resolution";
0292 mehCSCResRDPhi = iBooker.book1D(hcharname, hchartitle, 200, -0.2, 0.2);
0293 hcharname = "hRPCResX";
0294 hchartitle = "RPC rechits x resolution";
0295 mehRPCResX = iBooker.book1D(hcharname, hchartitle, 50, -5., 5.);
0296 }
0297
0298 void GlobalRecHitsAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0299 std::string MsgLoggerCat = "GlobalRecHitsAnalyzer_analyze";
0300
0301
0302 ++count;
0303
0304
0305 edm::RunNumber_t nrun = iEvent.id().run();
0306 edm::EventNumber_t nevt = iEvent.id().event();
0307
0308 if (verbosity > 0) {
0309 edm::LogInfo(MsgLoggerCat) << "Processing run " << nrun << ", event " << nevt << " (" << count << " events total)";
0310 } else if (verbosity == 0) {
0311 if (nevt % frequency == 0 || nevt == 1) {
0312 edm::LogInfo(MsgLoggerCat) << "Processing run " << nrun << ", event " << nevt << " (" << count
0313 << " events total)";
0314 }
0315 }
0316
0317
0318 if (getAllProvenances) {
0319 std::vector<const edm::StableProvenance*> AllProv;
0320 iEvent.getAllStableProvenance(AllProv);
0321
0322 if (verbosity >= 0)
0323 edm::LogInfo(MsgLoggerCat) << "Number of Provenances = " << AllProv.size();
0324
0325 if (printProvenanceInfo && (verbosity >= 0)) {
0326 TString eventout("\nProvenance info:\n");
0327
0328 for (unsigned int i = 0; i < AllProv.size(); ++i) {
0329 eventout += "\n ******************************";
0330 eventout += "\n Module : ";
0331 eventout += AllProv[i]->moduleLabel();
0332 eventout += "\n ProductID : ";
0333 eventout += AllProv[i]->productID().id();
0334 eventout += "\n ClassName : ";
0335 eventout += AllProv[i]->className();
0336 eventout += "\n InstanceName : ";
0337 eventout += AllProv[i]->productInstanceName();
0338 eventout += "\n BranchName : ";
0339 eventout += AllProv[i]->branchName();
0340 }
0341 eventout += "\n ******************************\n";
0342 edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0343 printProvenanceInfo = false;
0344 }
0345 getAllProvenances = false;
0346 }
0347
0348
0349
0350 fillECal(iEvent, iSetup);
0351
0352 fillHCal(iEvent, iSetup);
0353
0354 fillTrk(iEvent, iSetup);
0355
0356 fillMuon(iEvent, iSetup);
0357
0358 if (verbosity > 0)
0359 edm::LogInfo(MsgLoggerCat) << "Done gathering data from event.";
0360
0361 return;
0362 }
0363
0364 void GlobalRecHitsAnalyzer::fillECal(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0365 std::string MsgLoggerCat = "GlobalRecHitsAnalyzer_fillECal";
0366
0367 TString eventout;
0368 if (verbosity > 0)
0369 eventout = "\nGathering info:";
0370
0371
0372 edm::Handle<CrossingFrame<PCaloHit>> crossingFrame;
0373
0374
0375
0376
0377 edm::Handle<EBUncalibratedRecHitCollection> EcalUncalibRecHitEB;
0378 iEvent.getByToken(ECalUncalEBSrc_Token_, EcalUncalibRecHitEB);
0379 bool validUncalibRecHitEB = EcalUncalibRecHitEB.isValid();
0380 if (!validUncalibRecHitEB) {
0381 LogDebug(MsgLoggerCat) << "Unable to find EcalUncalRecHitEB in event!";
0382 }
0383
0384 edm::Handle<EBRecHitCollection> EcalRecHitEB;
0385 iEvent.getByToken(ECalEBSrc_Token_, EcalRecHitEB);
0386 bool validRecHitEB = EcalRecHitEB.isValid();
0387 if (!validRecHitEB) {
0388 LogDebug(MsgLoggerCat) << "Unable to find EcalRecHitEB in event!";
0389 }
0390
0391
0392 iEvent.getByToken(EBHits_Token_, crossingFrame);
0393 bool validXFrame = crossingFrame.isValid();
0394 if (!validXFrame) {
0395 LogDebug(MsgLoggerCat) << "Unable to find cal barrel crossingFrame in event!";
0396 }
0397
0398 MapType ebSimMap;
0399 if (validXFrame) {
0400 const MixCollection<PCaloHit> barrelHits(crossingFrame.product());
0401
0402 for (auto const& iHit : barrelHits) {
0403 EBDetId ebid = EBDetId(iHit.id());
0404
0405 uint32_t crystid = ebid.rawId();
0406 ebSimMap[crystid] += iHit.energy();
0407 }
0408 }
0409
0410 int nEBRecHits = 0;
0411
0412 if (validUncalibRecHitEB && validRecHitEB) {
0413 const EBUncalibratedRecHitCollection* EBUncalibRecHit = EcalUncalibRecHitEB.product();
0414 const EBRecHitCollection* EBRecHit = EcalRecHitEB.product();
0415
0416 for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EBUncalibRecHit->begin();
0417 uncalibRecHit != EBUncalibRecHit->end();
0418 ++uncalibRecHit) {
0419 EBDetId EBid = EBDetId(uncalibRecHit->id());
0420
0421 EcalRecHitCollection::const_iterator myRecHit = EBRecHit->find(EBid);
0422
0423 if (myRecHit != EBRecHit->end()) {
0424 ++nEBRecHits;
0425 mehEcalRes[1]->Fill(myRecHit->energy() - ebSimMap[EBid.rawId()]);
0426 }
0427 }
0428
0429 if (verbosity > 1) {
0430 eventout += "\n Number of EBRecHits collected:............ ";
0431 eventout += nEBRecHits;
0432 }
0433 mehEcaln[1]->Fill((float)nEBRecHits);
0434 }
0435
0436
0437
0438
0439 edm::Handle<EEUncalibratedRecHitCollection> EcalUncalibRecHitEE;
0440 iEvent.getByToken(ECalUncalEESrc_Token_, EcalUncalibRecHitEE);
0441 bool validuncalibRecHitEE = EcalUncalibRecHitEE.isValid();
0442 if (!validuncalibRecHitEE) {
0443 LogDebug(MsgLoggerCat) << "Unable to find EcalUncalRecHitEE in event!";
0444 }
0445
0446 edm::Handle<EERecHitCollection> EcalRecHitEE;
0447 iEvent.getByToken(ECalEESrc_Token_, EcalRecHitEE);
0448 bool validRecHitEE = EcalRecHitEE.isValid();
0449 if (!validRecHitEE) {
0450 LogDebug(MsgLoggerCat) << "Unable to find EcalRecHitEE in event!";
0451 }
0452
0453
0454 iEvent.getByToken(EEHits_Token_, crossingFrame);
0455 validXFrame = crossingFrame.isValid();
0456 if (!validXFrame) {
0457 LogDebug(MsgLoggerCat) << "Unable to find cal endcap crossingFrame in event!";
0458 }
0459
0460 MapType eeSimMap;
0461 if (validXFrame) {
0462 const MixCollection<PCaloHit> endcapHits(crossingFrame.product());
0463
0464 for (auto const& iHit : endcapHits) {
0465 EEDetId eeid = EEDetId(iHit.id());
0466
0467 uint32_t crystid = eeid.rawId();
0468 eeSimMap[crystid] += iHit.energy();
0469 }
0470 }
0471
0472 int nEERecHits = 0;
0473 if (validuncalibRecHitEE && validRecHitEE) {
0474
0475 const EEUncalibratedRecHitCollection* EEUncalibRecHit = EcalUncalibRecHitEE.product();
0476 const EERecHitCollection* EERecHit = EcalRecHitEE.product();
0477
0478 for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EEUncalibRecHit->begin();
0479 uncalibRecHit != EEUncalibRecHit->end();
0480 ++uncalibRecHit) {
0481 EEDetId EEid = EEDetId(uncalibRecHit->id());
0482
0483 EcalRecHitCollection::const_iterator myRecHit = EERecHit->find(EEid);
0484
0485 if (myRecHit != EERecHit->end()) {
0486 ++nEERecHits;
0487 mehEcalRes[0]->Fill(myRecHit->energy() - eeSimMap[EEid.rawId()]);
0488 }
0489 }
0490
0491 if (verbosity > 1) {
0492 eventout += "\n Number of EERecHits collected:............ ";
0493 eventout += nEERecHits;
0494 }
0495 mehEcaln[0]->Fill((float)nEERecHits);
0496 }
0497
0498
0499
0500
0501 edm::Handle<ESRecHitCollection> EcalRecHitES;
0502 iEvent.getByToken(ECalESSrc_Token_, EcalRecHitES);
0503 bool validRecHitES = EcalRecHitES.isValid();
0504 if (!validRecHitES) {
0505 LogDebug(MsgLoggerCat) << "Unable to find EcalRecHitES in event!";
0506 }
0507
0508
0509 iEvent.getByToken(ESHits_Token_, crossingFrame);
0510 validXFrame = crossingFrame.isValid();
0511 if (!validXFrame) {
0512 LogDebug(MsgLoggerCat) << "Unable to find cal preshower crossingFrame in event!";
0513 }
0514
0515 MapType esSimMap;
0516 if (validXFrame) {
0517 const MixCollection<PCaloHit> preshowerHits(crossingFrame.product());
0518
0519 for (auto const& iHit : preshowerHits) {
0520 ESDetId esid = ESDetId(iHit.id());
0521
0522 uint32_t crystid = esid.rawId();
0523 esSimMap[crystid] += iHit.energy();
0524 }
0525 }
0526
0527 int nESRecHits = 0;
0528 if (validRecHitES) {
0529
0530 const ESRecHitCollection* ESRecHit = EcalRecHitES.product();
0531 for (EcalRecHitCollection::const_iterator recHit = ESRecHit->begin(); recHit != ESRecHit->end(); ++recHit) {
0532 ESDetId ESid = ESDetId(recHit->id());
0533
0534 ++nESRecHits;
0535 mehEcalRes[2]->Fill(recHit->energy() - esSimMap[ESid.rawId()]);
0536 }
0537
0538 if (verbosity > 1) {
0539 eventout += "\n Number of ESRecHits collected:............ ";
0540 eventout += nESRecHits;
0541 }
0542 mehEcaln[2]->Fill(float(nESRecHits));
0543 }
0544
0545 if (verbosity > 0)
0546 edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0547
0548 return;
0549 }
0550
0551 void GlobalRecHitsAnalyzer::fillHCal(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0552 std::string MsgLoggerCat = "GlobalRecHitsAnalyzer_fillHCal";
0553
0554 TString eventout;
0555 if (verbosity > 0)
0556 eventout = "\nGathering info:";
0557
0558
0559 const auto& geometry = iSetup.getHandle(caloGeomToken_);
0560 if (!geometry.isValid()) {
0561 edm::LogWarning(MsgLoggerCat) << "Unable to find CaloGeometry in event!";
0562 return;
0563 }
0564
0565
0566
0567
0568 edm::Handle<edm::PCaloHitContainer> hcalHits;
0569 iEvent.getByToken(HCalSrc_Token_, hcalHits);
0570 bool validhcalHits = hcalHits.isValid();
0571 if (!validhcalHits) {
0572 LogDebug(MsgLoggerCat) << "Unable to find hcalHits in event!";
0573 }
0574
0575 std::map<HcalDetId, float> fHBEnergySimHits;
0576 std::map<HcalDetId, float> fHEEnergySimHits;
0577 std::map<HcalDetId, float> fHOEnergySimHits;
0578 std::map<HcalDetId, float> fHFEnergySimHits;
0579 if (validhcalHits) {
0580 const edm::PCaloHitContainer* simhitResult = hcalHits.product();
0581
0582 for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin(); simhits != simhitResult->end();
0583 ++simhits) {
0584 HcalDetId detId(simhits->id());
0585
0586 if (detId.subdet() == sdHcalBrl) {
0587 fHBEnergySimHits[detId] += simhits->energy();
0588 }
0589 if (detId.subdet() == sdHcalEC) {
0590 fHEEnergySimHits[detId] += simhits->energy();
0591 }
0592 if (detId.subdet() == sdHcalOut) {
0593 fHOEnergySimHits[detId] += simhits->energy();
0594 }
0595 if (detId.subdet() == sdHcalFwd) {
0596 fHFEnergySimHits[detId] += simhits->energy();
0597 }
0598 }
0599 }
0600
0601
0602
0603
0604 std::vector<edm::Handle<HBHERecHitCollection>> hbhe;
0605 iEvent.getManyByType(hbhe);
0606
0607 bool validHBHE = hbhe[0].isValid();
0608
0609 if (!validHBHE) {
0610 LogDebug(MsgLoggerCat) << "Unable to find any HBHERecHitCollections in event!";
0611 }
0612
0613 else {
0614 std::vector<edm::Handle<HBHERecHitCollection>>::iterator ihbhe;
0615 int iHB = 0;
0616 int iHE = 0;
0617 for (ihbhe = hbhe.begin(); ihbhe != hbhe.end(); ++ihbhe) {
0618 for (HBHERecHitCollection::const_iterator jhbhe = (*ihbhe)->begin(); jhbhe != (*ihbhe)->end(); ++jhbhe) {
0619 HcalDetId cell(jhbhe->id());
0620
0621 if (cell.subdet() == sdHcalBrl) {
0622 ++iHB;
0623 mehHcalRes[0]->Fill(jhbhe->energy() - fHBEnergySimHits[cell]);
0624 }
0625
0626 else if (cell.subdet() == sdHcalEC) {
0627 ++iHE;
0628 mehHcalRes[1]->Fill(jhbhe->energy() - fHEEnergySimHits[cell]);
0629 }
0630 }
0631 }
0632
0633 if (verbosity > 1) {
0634 eventout += "\n Number of HBRecHits collected:............ ";
0635 eventout += iHB;
0636 }
0637
0638 if (verbosity > 1) {
0639 eventout += "\n Number of HERecHits collected:............ ";
0640 eventout += iHE;
0641 }
0642 mehHcaln[0]->Fill((float)iHB);
0643 mehHcaln[1]->Fill((float)iHE);
0644 }
0645
0646
0647
0648
0649 std::vector<edm::Handle<HFRecHitCollection>> hf;
0650 iEvent.getManyByType(hf);
0651 bool validHF = hf[0].isValid();
0652 if (!validHF) {
0653 LogDebug(MsgLoggerCat) << "Unable to find any HFRecHitCollections in event!";
0654 } else {
0655 std::vector<edm::Handle<HFRecHitCollection>>::iterator ihf;
0656
0657 int iHF = 0;
0658 for (ihf = hf.begin(); ihf != hf.end(); ++ihf) {
0659 for (HFRecHitCollection::const_iterator jhf = (*ihf)->begin(); jhf != (*ihf)->end(); ++jhf) {
0660 HcalDetId cell(jhf->id());
0661
0662 if (cell.subdet() == sdHcalFwd) {
0663 ++iHF;
0664 mehHcalRes[2]->Fill(jhf->energy() - fHFEnergySimHits[cell]);
0665 }
0666 }
0667 }
0668
0669 if (verbosity > 1) {
0670 eventout += "\n Number of HFDigis collected:.............. ";
0671 eventout += iHF;
0672 }
0673 mehHcaln[2]->Fill((float)iHF);
0674 }
0675
0676
0677
0678
0679 std::vector<edm::Handle<HORecHitCollection>> ho;
0680 iEvent.getManyByType(ho);
0681 bool validHO = ho[0].isValid();
0682 if (!validHO) {
0683 LogDebug(MsgLoggerCat) << "Unable to find any HORecHitCollections in event!";
0684 } else {
0685 std::vector<edm::Handle<HORecHitCollection>>::iterator iho;
0686
0687 int iHO = 0;
0688 for (iho = ho.begin(); iho != ho.end(); ++iho) {
0689 for (HORecHitCollection::const_iterator jho = (*iho)->begin(); jho != (*iho)->end(); ++jho) {
0690 HcalDetId cell(jho->id());
0691
0692 if (cell.subdet() == sdHcalOut) {
0693 ++iHO;
0694 mehHcalRes[3]->Fill(jho->energy() - fHOEnergySimHits[cell]);
0695 }
0696 }
0697 }
0698
0699 if (verbosity > 1) {
0700 eventout += "\n Number of HODigis collected:.............. ";
0701 eventout += iHO;
0702 }
0703 mehHcaln[3]->Fill((float)iHO);
0704 }
0705
0706 if (verbosity > 0)
0707 edm::LogInfo(MsgLoggerCat) << eventout << "\n";
0708
0709 return;
0710 }
0711
0712 void GlobalRecHitsAnalyzer::fillTrk(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0713
0714 const TrackerTopology* const tTopo = &iSetup.getData(tTopoToken_);
0715 std::string MsgLoggerCat = "GlobalRecHitsAnalyzer_fillTrk";
0716 TString eventout;
0717 if (verbosity > 0)
0718 eventout = "\nGathering info:";
0719
0720
0721 edm::Handle<SiStripMatchedRecHit2DCollection> rechitsmatched;
0722 iEvent.getByToken(SiStripSrc_Token_, rechitsmatched);
0723 bool validstrip = rechitsmatched.isValid();
0724 if (!validstrip) {
0725 LogDebug(MsgLoggerCat) << "Unable to find stripmatchedrechits in event!";
0726 }
0727
0728 TrackerHitAssociator associate(iEvent, trackerHitAssociatorConfig_);
0729
0730 const auto& tGeomHandle = iSetup.getHandle(tGeomToken_);
0731 if (!tGeomHandle.isValid()) {
0732 edm::LogWarning(MsgLoggerCat) << "Unable to find TrackerDigiGeometry in event!";
0733 return;
0734 }
0735 const TrackerGeometry& tracker(*tGeomHandle);
0736
0737 if (validstrip) {
0738 int nStripBrl = 0, nStripFwd = 0;
0739
0740
0741 for (TrackerGeometry::DetContainer::const_iterator it = tGeomHandle->dets().begin();
0742 it != tGeomHandle->dets().end();
0743 ++it) {
0744 uint32_t myid = ((*it)->geographicalId()).rawId();
0745 DetId detid = ((*it)->geographicalId());
0746
0747
0748 SiStripMatchedRecHit2DCollection::const_iterator rechitmatchedMatch = rechitsmatched->find(detid);
0749
0750 if (rechitmatchedMatch != rechitsmatched->end()) {
0751 SiStripMatchedRecHit2DCollection::DetSet rechitmatchedRange = *rechitmatchedMatch;
0752 SiStripMatchedRecHit2DCollection::DetSet::const_iterator rechitmatchedRangeIteratorBegin =
0753 rechitmatchedRange.begin();
0754 SiStripMatchedRecHit2DCollection::DetSet::const_iterator rechitmatchedRangeIteratorEnd =
0755 rechitmatchedRange.end();
0756 SiStripMatchedRecHit2DCollection::DetSet::const_iterator itermatched;
0757
0758 for (itermatched = rechitmatchedRangeIteratorBegin; itermatched != rechitmatchedRangeIteratorEnd;
0759 ++itermatched) {
0760 SiStripMatchedRecHit2D const rechit = *itermatched;
0761 LocalPoint position = rechit.localPosition();
0762
0763 float mindist = 999999.;
0764 float distx = 999999.;
0765 float disty = 999999.;
0766 float dist = 999999.;
0767 std::pair<LocalPoint, LocalVector> closestPair;
0768 matched.clear();
0769
0770 float rechitmatchedx = position.x();
0771 float rechitmatchedy = position.y();
0772
0773 matched = associate.associateHit(rechit);
0774
0775 if (!matched.empty()) {
0776
0777 const GluedGeomDet* gluedDet = (const GluedGeomDet*)tracker.idToDet(rechit.geographicalId());
0778 const StripGeomDetUnit* partnerstripdet = (StripGeomDetUnit*)gluedDet->stereoDet();
0779 std::pair<LocalPoint, LocalVector> hitPair;
0780
0781 for (std::vector<PSimHit>::const_iterator m = matched.begin(); m != matched.end(); m++) {
0782
0783 hitPair = projectHit((*m), partnerstripdet, gluedDet->surface());
0784 distx = fabs(rechitmatchedx - hitPair.first.x());
0785 disty = fabs(rechitmatchedy - hitPair.first.y());
0786 dist = sqrt(distx * distx + disty * disty);
0787
0788 if (dist < mindist) {
0789 mindist = dist;
0790 closestPair = hitPair;
0791 }
0792 }
0793
0794
0795 if (detid.subdetId() == sdSiTIB) {
0796 ++nStripBrl;
0797
0798 if (tTopo->tibLayer(myid) == 1) {
0799 mehSiStripResX[8]->Fill(rechitmatchedx - closestPair.first.x());
0800 mehSiStripResY[8]->Fill(rechitmatchedy - closestPair.first.y());
0801 }
0802 if (tTopo->tibLayer(myid) == 2) {
0803 mehSiStripResX[9]->Fill(rechitmatchedx - closestPair.first.x());
0804 mehSiStripResY[9]->Fill(rechitmatchedy - closestPair.first.y());
0805 }
0806 if (tTopo->tibLayer(myid) == 3) {
0807 mehSiStripResX[10]->Fill(rechitmatchedx - closestPair.first.x());
0808 mehSiStripResY[10]->Fill(rechitmatchedy - closestPair.first.y());
0809 }
0810 if (tTopo->tibLayer(myid) == 4) {
0811 mehSiStripResX[11]->Fill(rechitmatchedx - closestPair.first.x());
0812 mehSiStripResY[11]->Fill(rechitmatchedy - closestPair.first.y());
0813 }
0814 }
0815
0816
0817 if (detid.subdetId() == sdSiTOB) {
0818 ++nStripBrl;
0819
0820 if (tTopo->tobLayer(myid) == 1) {
0821 mehSiStripResX[15]->Fill(rechitmatchedx - closestPair.first.x());
0822 mehSiStripResY[15]->Fill(rechitmatchedy - closestPair.first.y());
0823 }
0824 if (tTopo->tobLayer(myid) == 2) {
0825 mehSiStripResX[16]->Fill(rechitmatchedx - closestPair.first.x());
0826 mehSiStripResY[16]->Fill(rechitmatchedy - closestPair.first.y());
0827 }
0828 if (tTopo->tobLayer(myid) == 3) {
0829 mehSiStripResX[17]->Fill(rechitmatchedx - closestPair.first.x());
0830 mehSiStripResY[17]->Fill(rechitmatchedy - closestPair.first.y());
0831 }
0832 if (tTopo->tobLayer(myid) == 4) {
0833 mehSiStripResX[18]->Fill(rechitmatchedx - closestPair.first.x());
0834 mehSiStripResY[18]->Fill(rechitmatchedy - closestPair.first.y());
0835 }
0836 }
0837
0838
0839 if (detid.subdetId() == sdSiTID) {
0840 ++nStripFwd;
0841
0842 if (tTopo->tidWheel(myid) == 1) {
0843 mehSiStripResX[12]->Fill(rechitmatchedx - closestPair.first.x());
0844 mehSiStripResY[12]->Fill(rechitmatchedy - closestPair.first.y());
0845 }
0846 if (tTopo->tidWheel(myid) == 2) {
0847 mehSiStripResX[13]->Fill(rechitmatchedx - closestPair.first.x());
0848 mehSiStripResY[13]->Fill(rechitmatchedy - closestPair.first.y());
0849 }
0850 if (tTopo->tidWheel(myid) == 3) {
0851 mehSiStripResX[14]->Fill(rechitmatchedx - closestPair.first.x());
0852 mehSiStripResY[14]->Fill(rechitmatchedy - closestPair.first.y());
0853 }
0854 }
0855
0856
0857 if (detid.subdetId() == sdSiTEC) {
0858 ++nStripFwd;
0859
0860 if (tTopo->tecWheel(myid) == 1) {
0861 mehSiStripResX[0]->Fill(rechitmatchedx - closestPair.first.x());
0862 mehSiStripResY[0]->Fill(rechitmatchedy - closestPair.first.y());
0863 }
0864 if (tTopo->tecWheel(myid) == 2) {
0865 mehSiStripResX[1]->Fill(rechitmatchedx - closestPair.first.x());
0866 mehSiStripResY[1]->Fill(rechitmatchedy - closestPair.first.y());
0867 }
0868 if (tTopo->tecWheel(myid) == 3) {
0869 mehSiStripResX[2]->Fill(rechitmatchedx - closestPair.first.x());
0870 mehSiStripResY[2]->Fill(rechitmatchedy - closestPair.first.y());
0871 }
0872 if (tTopo->tecWheel(myid) == 4) {
0873 mehSiStripResX[3]->Fill(rechitmatchedx - closestPair.first.x());
0874 mehSiStripResY[3]->Fill(rechitmatchedy - closestPair.first.y());
0875 }
0876 if (tTopo->tecWheel(myid) == 5) {
0877 mehSiStripResX[4]->Fill(rechitmatchedx - closestPair.first.x());
0878 mehSiStripResY[4]->Fill(rechitmatchedy - closestPair.first.y());
0879 }
0880 if (tTopo->tecWheel(myid) == 6) {
0881 mehSiStripResX[5]->Fill(rechitmatchedx - closestPair.first.x());
0882 mehSiStripResY[5]->Fill(rechitmatchedy - closestPair.first.y());
0883 }
0884 if (tTopo->tecWheel(myid) == 7) {
0885 mehSiStripResX[6]->Fill(rechitmatchedx - closestPair.first.x());
0886 mehSiStripResY[6]->Fill(rechitmatchedy - closestPair.first.y());
0887 }
0888 if (tTopo->tecWheel(myid) == 8) {
0889 mehSiStripResX[7]->Fill(rechitmatchedx - closestPair.first.x());
0890 mehSiStripResY[7]->Fill(rechitmatchedy - closestPair.first.y());
0891 }
0892 }
0893
0894 }
0895 }
0896 }
0897 }
0898
0899 if (verbosity > 1) {
0900 eventout += "\n Number of BrlStripRecHits collected:...... ";
0901 eventout += nStripBrl;
0902 }
0903
0904 for (int i = 8; i < 12; ++i) {
0905 mehSiStripn[i]->Fill((float)nStripBrl);
0906 }
0907 for (int i = 16; i < 19; ++i) {
0908 mehSiStripn[i]->Fill((float)nStripBrl);
0909 }
0910
0911 if (verbosity > 1) {
0912 eventout += "\n Number of FrwdStripRecHits collected:..... ";
0913 eventout += nStripFwd;
0914 }
0915 for (int i = 0; i < 8; ++i) {
0916 mehSiStripn[i]->Fill((float)nStripFwd);
0917 }
0918 for (int i = 12; i < 16; ++i) {
0919 mehSiStripn[i]->Fill((float)nStripFwd);
0920 }
0921 }
0922
0923
0924
0925 edm::Handle<SiPixelRecHitCollection> recHitColl;
0926 iEvent.getByToken(SiPxlSrc_Token_, recHitColl);
0927 bool validpixel = recHitColl.isValid();
0928 if (!validpixel) {
0929 LogDebug(MsgLoggerCat) << "Unable to find SiPixelRecHitCollection in event!";
0930 } else {
0931 int nPxlBrl = 0, nPxlFwd = 0;
0932
0933 for (TrackerGeometry::DetContainer::const_iterator it = tGeomHandle->dets().begin();
0934 it != tGeomHandle->dets().end();
0935 ++it) {
0936 uint32_t myid = ((*it)->geographicalId()).rawId();
0937 DetId detId = ((*it)->geographicalId());
0938 int subid = detId.subdetId();
0939
0940 if (!((subid == sdPxlBrl) || (subid == sdPxlFwd)))
0941 continue;
0942
0943 SiPixelRecHitCollection::const_iterator pixeldet = recHitColl->find(detId);
0944 if (pixeldet == recHitColl->end())
0945 continue;
0946 SiPixelRecHitCollection::DetSet pixelrechitRange = *pixeldet;
0947 SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorBegin = pixelrechitRange.begin();
0948 SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorEnd = pixelrechitRange.end();
0949 SiPixelRecHitCollection::DetSet::const_iterator pixeliter = pixelrechitRangeIteratorBegin;
0950
0951 std::vector<PSimHit> matched;
0952
0953
0954 for (; pixeliter != pixelrechitRangeIteratorEnd; ++pixeliter) {
0955 matched.clear();
0956 matched = associate.associateHit(*pixeliter);
0957
0958 if (!matched.empty()) {
0959 float closest = 9999.9;
0960 LocalPoint lp = pixeliter->localPosition();
0961 float rechit_x = lp.x();
0962 float rechit_y = lp.y();
0963
0964 float sim_x = 0.;
0965 float sim_y = 0.;
0966
0967
0968 for (std::vector<PSimHit>::const_iterator m = matched.begin(); m != matched.end(); ++m) {
0969 float sim_x1 = (*m).entryPoint().x();
0970 float sim_x2 = (*m).exitPoint().x();
0971 float sim_xpos = 0.5 * (sim_x1 + sim_x2);
0972
0973 float sim_y1 = (*m).entryPoint().y();
0974 float sim_y2 = (*m).exitPoint().y();
0975 float sim_ypos = 0.5 * (sim_y1 + sim_y2);
0976
0977 float x_res = fabs(sim_xpos - rechit_x);
0978 float y_res = fabs(sim_ypos - rechit_y);
0979
0980 float dist = sqrt(x_res * x_res + y_res * y_res);
0981
0982 if (dist < closest) {
0983 closest = dist;
0984 sim_x = sim_xpos;
0985 sim_y = sim_ypos;
0986 }
0987 }
0988
0989
0990 if (subid == sdPxlBrl) {
0991 ++nPxlBrl;
0992
0993 if (tTopo->pxbLayer(myid) == 1) {
0994 mehSiPixelResX[0]->Fill(rechit_x - sim_x);
0995 mehSiPixelResY[0]->Fill(rechit_y - sim_y);
0996 }
0997 if (tTopo->pxbLayer(myid) == 2) {
0998 mehSiPixelResX[1]->Fill(rechit_x - sim_x);
0999 mehSiPixelResY[1]->Fill(rechit_y - sim_y);
1000 }
1001 if (tTopo->pxbLayer(myid) == 3) {
1002 mehSiPixelResX[2]->Fill(rechit_x - sim_x);
1003 mehSiPixelResY[2]->Fill(rechit_y - sim_y);
1004 }
1005 }
1006
1007
1008 if (subid == sdPxlFwd) {
1009 ++nPxlFwd;
1010
1011 if (tTopo->pxfDisk(myid) == 1) {
1012 if (tTopo->pxfSide(myid) == 1) {
1013 mehSiPixelResX[3]->Fill(rechit_x - sim_x);
1014 mehSiPixelResY[3]->Fill(rechit_y - sim_y);
1015 }
1016 if (tTopo->pxfSide(myid) == 2) {
1017 mehSiPixelResX[4]->Fill(rechit_x - sim_x);
1018 mehSiPixelResY[4]->Fill(rechit_y - sim_y);
1019 }
1020 }
1021 if (tTopo->pxfDisk(myid) == 2) {
1022 if (tTopo->pxfSide(myid) == 1) {
1023 mehSiPixelResX[5]->Fill(rechit_x - sim_x);
1024 mehSiPixelResY[5]->Fill(rechit_y - sim_y);
1025 }
1026 if (tTopo->pxfSide(myid) == 2) {
1027 mehSiPixelResX[6]->Fill(rechit_x - sim_x);
1028 mehSiPixelResY[6]->Fill(rechit_y - sim_y);
1029 }
1030 }
1031 }
1032 }
1033 }
1034 }
1035
1036 if (verbosity > 1) {
1037 eventout += "\n Number of BrlPixelRecHits collected:...... ";
1038 eventout += nPxlBrl;
1039 }
1040 for (int i = 0; i < 3; ++i) {
1041 mehSiPixeln[i]->Fill((float)nPxlBrl);
1042 }
1043
1044 if (verbosity > 1) {
1045 eventout += "\n Number of FrwdPixelRecHits collected:..... ";
1046 eventout += nPxlFwd;
1047 }
1048
1049 for (int i = 3; i < 7; ++i) {
1050 mehSiPixeln[i]->Fill((float)nPxlFwd);
1051 }
1052 }
1053
1054 if (verbosity > 0)
1055 edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1056
1057 return;
1058 }
1059
1060 void GlobalRecHitsAnalyzer::fillMuon(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
1061 std::string MsgLoggerCat = "GlobalRecHitsAnalyzer_fillMuon";
1062
1063 TString eventout;
1064 if (verbosity > 0)
1065 eventout = "\nGathering info:";
1066
1067
1068 const auto& dtGeom = iSetup.getHandle(dtGeomToken_);
1069 if (!dtGeom.isValid()) {
1070 edm::LogWarning(MsgLoggerCat) << "Unable to find DTMuonGeometryRecord in event!";
1071 return;
1072 }
1073
1074 edm::Handle<edm::PSimHitContainer> dtsimHits;
1075 iEvent.getByToken(MuDTSimSrc_Token_, dtsimHits);
1076 bool validdtsim = dtsimHits.isValid();
1077 if (!validdtsim) {
1078 LogDebug(MsgLoggerCat) << "Unable to find dtsimHits in event!";
1079 }
1080
1081 edm::Handle<DTRecHitCollection> dtRecHits;
1082 iEvent.getByToken(MuDTSrc_Token_, dtRecHits);
1083 bool validdtrec = dtRecHits.isValid();
1084 if (!validdtrec) {
1085 LogDebug(MsgLoggerCat) << "Unable to find dtRecHits in event!";
1086 }
1087
1088 if (validdtsim && validdtrec) {
1089 std::map<DTWireId, edm::PSimHitContainer> simHitsPerWire =
1090 DTHitQualityUtils::mapSimHitsPerWire(*(dtsimHits.product()));
1091
1092 std::map<DTWireId, std::vector<DTRecHit1DPair>> recHitsPerWire = map1DRecHitsPerWire(dtRecHits.product());
1093
1094 int nDt = compute(dtGeom.product(), simHitsPerWire, recHitsPerWire, 1);
1095
1096 if (verbosity > 1) {
1097 eventout += "\n Number of DtMuonRecHits collected:........ ";
1098 eventout += nDt;
1099 }
1100 mehDtMuonn->Fill(float(nDt));
1101 }
1102
1103
1104
1105 theMap.clear();
1106 edm::Handle<CrossingFrame<PSimHit>> cf;
1107
1108 iEvent.getByToken(MuCSCHits_Token_, cf);
1109 bool validXFrame = cf.isValid();
1110 if (!validXFrame) {
1111 LogDebug(MsgLoggerCat) << "Unable to find muo CSC crossingFrame in event!";
1112 } else {
1113 const MixCollection<PSimHit> simHits(cf.product());
1114
1115
1116 for (auto const& iHit : simHits) {
1117 theMap[iHit.detUnitId()].push_back(iHit);
1118 }
1119 }
1120
1121
1122 const auto& hGeom = iSetup.getHandle(cscGeomToken_);
1123 if (!hGeom.isValid()) {
1124 edm::LogWarning(MsgLoggerCat) << "Unable to find CSCMuonGeometryRecord in event!";
1125 return;
1126 }
1127 const CSCGeometry* theCSCGeometry = &*hGeom;
1128
1129
1130 edm::Handle<CSCRecHit2DCollection> hRecHits;
1131 iEvent.getByToken(MuCSCSrc_Token_, hRecHits);
1132 bool validCSC = hRecHits.isValid();
1133 if (!validCSC) {
1134 LogDebug(MsgLoggerCat) << "Unable to find CSC RecHits in event!";
1135 } else {
1136 const CSCRecHit2DCollection* cscRecHits = hRecHits.product();
1137
1138 int nCSC = 0;
1139 for (CSCRecHit2DCollection::const_iterator recHitItr = cscRecHits->begin(); recHitItr != cscRecHits->end();
1140 ++recHitItr) {
1141 int detId = (*recHitItr).cscDetId().rawId();
1142
1143 edm::PSimHitContainer simHits;
1144 std::map<int, edm::PSimHitContainer>::const_iterator mapItr = theMap.find(detId);
1145 if (mapItr != theMap.end()) {
1146 simHits = mapItr->second;
1147 }
1148
1149 if (simHits.size() == 1) {
1150 ++nCSC;
1151
1152 const GeomDetUnit* detUnit = theCSCGeometry->idToDetUnit(CSCDetId(detId));
1153 const CSCLayer* layer = dynamic_cast<const CSCLayer*>(detUnit);
1154
1155 int chamberType = layer->chamber()->specs()->chamberType();
1156 plotResolution(simHits[0], *recHitItr, layer, chamberType);
1157 }
1158 }
1159
1160 if (verbosity > 1) {
1161 eventout += "\n Number of CSCRecHits collected:........... ";
1162 eventout += nCSC;
1163 }
1164 mehCSCn->Fill((float)nCSC);
1165 }
1166
1167
1168 std::map<double, int> mapsim, maprec;
1169 std::map<int, double> nmapsim, nmaprec;
1170 const auto& rpcGeom = iSetup.getHandle(rpcGeomToken_);
1171 if (!rpcGeom.isValid()) {
1172 edm::LogWarning(MsgLoggerCat) << "Unable to find RPCMuonGeometryRecord in event!";
1173 return;
1174 }
1175
1176 edm::Handle<edm::PSimHitContainer> simHit;
1177 iEvent.getByToken(MuRPCSimSrc_Token_, simHit);
1178 bool validrpcsim = simHit.isValid();
1179 if (!validrpcsim) {
1180 LogDebug(MsgLoggerCat) << "Unable to find RPCSimHit in event!";
1181 }
1182
1183 edm::Handle<RPCRecHitCollection> recHit;
1184 iEvent.getByToken(MuRPCSrc_Token_, recHit);
1185 bool validrpc = recHit.isValid();
1186 if (!validrpc) {
1187 LogDebug(MsgLoggerCat) << "Unable to find RPCRecHit in event!";
1188 }
1189
1190 if (validrpc) {
1191 int nRPC = 0;
1192 RPCRecHitCollection::const_iterator recIt;
1193 int nrec = 0;
1194 for (recIt = recHit->begin(); recIt != recHit->end(); ++recIt) {
1195 RPCDetId Rid = (RPCDetId)(*recIt).rpcId();
1196 const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(Rid));
1197 if (roll->isForward()) {
1198 if (verbosity > 1) {
1199 eventout += "\n Number of RPCRecHits collected:........... ";
1200 eventout += nRPC;
1201 }
1202
1203 if (verbosity > 0)
1204 edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1205 return;
1206 }
1207 nrec = nrec + 1;
1208 LocalPoint rhitlocal = (*recIt).localPosition();
1209 double rhitlocalx = rhitlocal.x();
1210 maprec[rhitlocalx] = nrec;
1211 }
1212
1213 int i = 0;
1214 for (std::map<double, int>::iterator iter = maprec.begin(); iter != maprec.end(); ++iter) {
1215 i = i + 1;
1216 nmaprec[i] = (*iter).first;
1217 }
1218
1219 int nsim = 0;
1220 if (validrpcsim) {
1221 edm::PSimHitContainer::const_iterator simIt;
1222 for (simIt = simHit->begin(); simIt != simHit->end(); simIt++) {
1223 int ptype = (*simIt).particleType();
1224 if (ptype == 13 || ptype == -13) {
1225 nsim = nsim + 1;
1226 LocalPoint shitlocal = (*simIt).localPosition();
1227 double shitlocalx = shitlocal.x();
1228 mapsim[shitlocalx] = nsim;
1229 }
1230 }
1231
1232 i = 0;
1233 for (std::map<double, int>::iterator iter = mapsim.begin(); iter != mapsim.end(); ++iter) {
1234 i = i + 1;
1235 nmapsim[i] = (*iter).first;
1236 }
1237 }
1238
1239 if (nsim == nrec) {
1240 for (int r = 0; r < nsim; r++) {
1241 ++nRPC;
1242 mehRPCResX->Fill(nmaprec[r + 1] - nmapsim[r + 1]);
1243 }
1244 }
1245
1246 if (verbosity > 1) {
1247 eventout += "\n Number of RPCRecHits collected:........... ";
1248 eventout += nRPC;
1249 }
1250 mehRPCn->Fill((float)nRPC);
1251 }
1252
1253 if (verbosity > 0)
1254 edm::LogInfo(MsgLoggerCat) << eventout << "\n";
1255
1256 return;
1257 }
1258
1259
1260 std::pair<LocalPoint, LocalVector> GlobalRecHitsAnalyzer::projectHit(const PSimHit& hit,
1261 const StripGeomDetUnit* stripDet,
1262 const BoundPlane& plane) {
1263 const StripTopology& topol = stripDet->specificTopology();
1264 GlobalPoint globalpos = stripDet->surface().toGlobal(hit.localPosition());
1265 LocalPoint localHit = plane.toLocal(globalpos);
1266
1267 LocalVector locdir = hit.localDirection();
1268
1269
1270 GlobalVector globaldir = stripDet->surface().toGlobal(locdir);
1271 LocalVector dir = plane.toLocal(globaldir);
1272 float scale = -localHit.z() / dir.z();
1273
1274 LocalPoint projectedPos = localHit + scale * dir;
1275
1276 float selfAngle = topol.stripAngle(topol.strip(hit.localPosition()));
1277
1278
1279 LocalVector stripDir(sin(selfAngle), cos(selfAngle), 0);
1280
1281 LocalVector localStripDir(plane.toLocal(stripDet->surface().toGlobal(stripDir)));
1282
1283 return std::pair<LocalPoint, LocalVector>(projectedPos, localStripDir);
1284 }
1285
1286
1287 std::map<DTWireId, std::vector<DTRecHit1DPair>> GlobalRecHitsAnalyzer::map1DRecHitsPerWire(
1288 const DTRecHitCollection* dt1DRecHitPairs) {
1289 std::map<DTWireId, std::vector<DTRecHit1DPair>> ret;
1290
1291 for (DTRecHitCollection::const_iterator rechit = dt1DRecHitPairs->begin(); rechit != dt1DRecHitPairs->end();
1292 rechit++) {
1293 ret[(*rechit).wireId()].push_back(*rechit);
1294 }
1295
1296 return ret;
1297 }
1298
1299
1300 float GlobalRecHitsAnalyzer::simHitDistFromWire(const DTLayer* layer, DTWireId wireId, const PSimHit& hit) {
1301 float xwire = layer->specificTopology().wirePosition(wireId.wire());
1302 LocalPoint entryP = hit.entryPoint();
1303 LocalPoint exitP = hit.exitPoint();
1304 float xEntry = entryP.x() - xwire;
1305 float xExit = exitP.x() - xwire;
1306
1307
1308 return fabs(xEntry - (entryP.z() * (xExit - xEntry)) / (exitP.z() - entryP.z()));
1309 }
1310
1311
1312 template <typename type>
1313 const type* GlobalRecHitsAnalyzer::findBestRecHit(const DTLayer* layer,
1314 DTWireId wireId,
1315 const std::vector<type>& recHits,
1316 const float simHitDist) {
1317 float res = 99999;
1318 const type* theBestRecHit = nullptr;
1319
1320 for (typename std::vector<type>::const_iterator recHit = recHits.begin(); recHit != recHits.end(); recHit++) {
1321 float distTmp = recHitDistFromWire(*recHit, layer);
1322 if (fabs(distTmp - simHitDist) < res) {
1323 res = fabs(distTmp - simHitDist);
1324 theBestRecHit = &(*recHit);
1325 }
1326 }
1327
1328 return theBestRecHit;
1329 }
1330
1331
1332 float GlobalRecHitsAnalyzer::recHitDistFromWire(const DTRecHit1DPair& hitPair, const DTLayer* layer) {
1333
1334 return fabs(hitPair.localPosition(DTEnums::Left).x() - hitPair.localPosition(DTEnums::Right).x()) / 2.;
1335 }
1336
1337
1338 float GlobalRecHitsAnalyzer::recHitDistFromWire(const DTRecHit1D& recHit, const DTLayer* layer) {
1339 return fabs(recHit.localPosition().x() - layer->specificTopology().wirePosition(recHit.wireId().wire()));
1340 }
1341
1342 template <typename type>
1343 int GlobalRecHitsAnalyzer::compute(const DTGeometry* dtGeom,
1344 const std::map<DTWireId, std::vector<PSimHit>>& _simHitsPerWire,
1345 const std::map<DTWireId, std::vector<type>>& _recHitsPerWire,
1346 int step) {
1347 std::map<DTWireId, std::vector<PSimHit>> simHitsPerWire = _simHitsPerWire;
1348 std::map<DTWireId, std::vector<type>> recHitsPerWire = _recHitsPerWire;
1349 int nDt = 0;
1350
1351 for (std::map<DTWireId, std::vector<PSimHit>>::const_iterator wireAndSHits = simHitsPerWire.begin();
1352 wireAndSHits != simHitsPerWire.end();
1353 wireAndSHits++) {
1354 DTWireId wireId = (*wireAndSHits).first;
1355 std::vector<PSimHit> simHitsInCell = (*wireAndSHits).second;
1356
1357
1358 const DTLayer* layer = dtGeom->layer(wireId);
1359
1360
1361 const PSimHit* muSimHit = DTHitQualityUtils::findMuSimHit(simHitsInCell);
1362 if (muSimHit == nullptr) {
1363 continue;
1364 }
1365
1366
1367 float simHitWireDist = simHitDistFromWire(layer, wireId, *muSimHit);
1368
1369 if (simHitWireDist > 2.1) {
1370 continue;
1371 }
1372
1373
1374 if (recHitsPerWire.find(wireId) == recHitsPerWire.end()) {
1375 continue;
1376 } else {
1377 std::vector<type> recHits = recHitsPerWire[wireId];
1378
1379
1380 const type* theBestRecHit = findBestRecHit(layer, wireId, recHits, simHitWireDist);
1381
1382 float recHitWireDist = recHitDistFromWire(*theBestRecHit, layer);
1383
1384 ++nDt;
1385
1386 mehDtMuonRes->Fill(recHitWireDist - simHitWireDist);
1387
1388 }
1389 }
1390
1391 return nDt;
1392 }
1393
1394 void GlobalRecHitsAnalyzer::plotResolution(const PSimHit& simHit,
1395 const CSCRecHit2D& recHit,
1396 const CSCLayer* layer,
1397 int chamberType) {
1398 GlobalPoint simHitPos = layer->toGlobal(simHit.localPosition());
1399 GlobalPoint recHitPos = layer->toGlobal(recHit.localPosition());
1400
1401 mehCSCResRDPhi->Fill(recHitPos.phi() - simHitPos.phi());
1402 }