21 getAllProvenances(
false),
22 printProvenanceInfo(
false),
23 trackerHitAssociatorConfig_(iPSet, consumesCollector()),
31 consumesMany<edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>>();
32 consumesMany<edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>>();
33 consumesMany<edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>>();
34 std::string MsgLoggerCat =
"GlobalRecHitsAnalyzer_GlobalRecHitsAnalyzer";
96 edm::LogInfo(MsgLoggerCat) <<
"\n===============================\n"
97 <<
"Initialized as EDProducer with parameter values:\n"
98 <<
" Name = " <<
fName <<
"\n"
122 <<
"===============================\n";
130 string SiStripString[19] = {
"TECW1",
149 for (
int i = 0;
i < 19; ++
i) {
154 string hcharname, hchartitle;
156 for (
int amend = 0; amend < 19; ++amend) {
157 hcharname =
"hSiStripn_" + SiStripString[amend];
158 hchartitle = SiStripString[amend] +
" rechits";
162 hcharname =
"hSiStripResX_" + SiStripString[amend];
163 hchartitle = SiStripString[amend] +
" rechit x resolution";
167 hcharname =
"hSiStripResY_" + SiStripString[amend];
168 hchartitle = SiStripString[amend] +
" rechit y resolution";
176 string HCalString[4] = {
"HB",
"HE",
"HF",
"HO"};
177 float HCalnUpper[4] = {3000., 3000., 3000., 3000.};
178 float HCalnLower[4] = {0., 0., 0., 0.};
179 for (
int j = 0;
j < 4; ++
j) {
185 for (
int amend = 0; amend < 4; ++amend) {
186 hcharname =
"hHcaln_" + HCalString[amend];
187 hchartitle = HCalString[amend] +
" rechits";
188 mehHcaln[amend] = iBooker.
book1D(hcharname, hchartitle, 1000, HCalnLower[amend], HCalnUpper[amend]);
191 hcharname =
"hHcalRes_" + HCalString[amend];
192 hchartitle = HCalString[amend] +
" rechit resolution";
199 string ECalString[3] = {
"EB",
"EE",
"ES"};
200 int ECalnBins[3] = {1000, 3000, 150};
201 float ECalnUpper[3] = {20000., 62000., 3000.};
202 float ECalnLower[3] = {0., 0., 0.};
203 int ECalResBins[3] = {200, 200, 200};
204 float ECalResUpper[3] = {1., 0.3, .0002};
205 float ECalResLower[3] = {-1., -0.3, -.0002};
206 for (
int i = 0;
i < 3; ++
i) {
212 for (
int amend = 0; amend < 3; ++amend) {
213 hcharname =
"hEcaln_" + ECalString[amend];
214 hchartitle = ECalString[amend] +
" rechits";
215 mehEcaln[amend] = iBooker.
book1D(hcharname, hchartitle, ECalnBins[amend], ECalnLower[amend], ECalnUpper[amend]);
218 hcharname =
"hEcalRes_" + ECalString[amend];
219 hchartitle = ECalString[amend] +
" rechit resolution";
221 iBooker.
book1D(hcharname, hchartitle, ECalResBins[amend], ECalResLower[amend], ECalResUpper[amend]);
227 string SiPixelString[7] = {
"BRL1",
"BRL2",
"BRL3",
"FWD1n",
"FWD1p",
"FWD2n",
"FWD2p"};
228 for (
int j = 0;
j < 7; ++
j) {
235 for (
int amend = 0; amend < 7; ++amend) {
236 hcharname =
"hSiPixeln_" + SiPixelString[amend];
237 hchartitle = SiPixelString[amend] +
" rechits";
241 hcharname =
"hSiPixelResX_" + SiPixelString[amend];
242 hchartitle = SiPixelString[amend] +
" rechit x resolution";
246 hcharname =
"hSiPixelResY_" + SiPixelString[amend];
247 hchartitle = SiPixelString[amend] +
" rechit y resolution";
261 string n_List[3] = {
"hDtMuonn",
"hCSCn",
"hRPCn"};
262 string hist_string[3] = {
"Dt",
"CSC",
"RPC"};
264 for (
int amend = 0; amend < 3; ++amend) {
265 hchartitle = hist_string[amend] +
" rechits";
272 mehCSCn = iBooker.
book1D(n_List[amend], hchartitle, 50, 0., 500.);
277 mehRPCn = iBooker.
book1D(n_List[amend], hchartitle, 50, 0., 500.);
287 hcharname =
"hDtMuonRes";
288 hchartitle =
"DT wire distance resolution";
290 hcharname =
"CSCResRDPhi";
291 hchartitle =
"CSC perp*dphi resolution";
293 hcharname =
"hRPCResX";
294 hchartitle =
"RPC rechits x resolution";
299 std::string MsgLoggerCat =
"GlobalRecHitsAnalyzer_analyze";
309 edm::LogInfo(MsgLoggerCat) <<
"Processing run " << nrun <<
", event " <<
nevt <<
" (" <<
count <<
" events total)";
319 std::vector<const edm::StableProvenance*> AllProv;
320 iEvent.getAllStableProvenance(AllProv);
323 edm::LogInfo(MsgLoggerCat) <<
"Number of Provenances = " << AllProv.size();
326 TString eventout(
"\nProvenance info:\n");
328 for (
unsigned int i = 0;
i < AllProv.size(); ++
i) {
329 eventout +=
"\n ******************************";
330 eventout +=
"\n Module : ";
331 eventout += AllProv[
i]->moduleLabel();
332 eventout +=
"\n ProductID : ";
333 eventout += AllProv[
i]->productID().id();
334 eventout +=
"\n ClassName : ";
335 eventout += AllProv[
i]->className();
336 eventout +=
"\n InstanceName : ";
337 eventout += AllProv[
i]->productInstanceName();
338 eventout +=
"\n BranchName : ";
339 eventout += AllProv[
i]->branchName();
341 eventout +=
"\n ******************************\n";
359 edm::LogInfo(MsgLoggerCat) <<
"Done gathering data from event.";
365 std::string MsgLoggerCat =
"GlobalRecHitsAnalyzer_fillECal";
369 eventout =
"\nGathering info:";
379 bool validUncalibRecHitEB =
true;
380 if (!EcalUncalibRecHitEB.
isValid()) {
381 LogDebug(MsgLoggerCat) <<
"Unable to find EcalUncalRecHitEB in event!";
382 validUncalibRecHitEB =
false;
387 bool validRecHitEB =
true;
389 LogDebug(MsgLoggerCat) <<
"Unable to find EcalRecHitEB in event!";
390 validRecHitEB =
false;
395 bool validXFrame =
true;
396 if (!crossingFrame.
isValid()) {
397 LogDebug(MsgLoggerCat) <<
"Unable to find cal barrel crossingFrame in event!";
408 uint32_t crystid = ebid.
rawId();
409 ebSimMap[crystid] += iHit.energy();
415 if (validUncalibRecHitEB && validRecHitEB) {
433 eventout +=
"\n Number of EBRecHits collected:............ ";
434 eventout += nEBRecHits;
444 bool validuncalibRecHitEE =
true;
445 if (!EcalUncalibRecHitEE.
isValid()) {
446 LogDebug(MsgLoggerCat) <<
"Unable to find EcalUncalRecHitEE in event!";
447 validuncalibRecHitEE =
false;
452 bool validRecHitEE =
true;
454 LogDebug(MsgLoggerCat) <<
"Unable to find EcalRecHitEE in event!";
455 validRecHitEE =
false;
461 if (!crossingFrame.
isValid()) {
462 LogDebug(MsgLoggerCat) <<
"Unable to find cal endcap crossingFrame in event!";
473 uint32_t crystid = eeid.
rawId();
474 eeSimMap[crystid] += iHit.energy();
479 if (validuncalibRecHitEE && validRecHitEE) {
498 eventout +=
"\n Number of EERecHits collected:............ ";
499 eventout += nEERecHits;
509 bool validRecHitES =
true;
511 LogDebug(MsgLoggerCat) <<
"Unable to find EcalRecHitES in event!";
512 validRecHitES =
false;
518 if (!crossingFrame.
isValid()) {
519 LogDebug(MsgLoggerCat) <<
"Unable to find cal preshower crossingFrame in event!";
530 uint32_t crystid = esid.
rawId();
531 esSimMap[crystid] += iHit.energy();
547 eventout +=
"\n Number of ESRecHits collected:............ ";
548 eventout += nESRecHits;
560 std::string MsgLoggerCat =
"GlobalRecHitsAnalyzer_fillHCal";
564 eventout =
"\nGathering info:";
569 edm::LogWarning(MsgLoggerCat) <<
"Unable to find CaloGeometry in event!";
578 bool validhcalHits =
true;
580 LogDebug(MsgLoggerCat) <<
"Unable to find hcalHits in event!";
581 validhcalHits =
false;
584 std::map<HcalDetId, float> fHBEnergySimHits;
585 std::map<HcalDetId, float> fHEEnergySimHits;
586 std::map<HcalDetId, float> fHOEnergySimHits;
587 std::map<HcalDetId, float> fHFEnergySimHits;
591 for (std::vector<PCaloHit>::const_iterator
simhits = simhitResult->begin();
simhits != simhitResult->end();
596 fHBEnergySimHits[detId] +=
simhits->energy();
599 fHEEnergySimHits[detId] +=
simhits->energy();
602 fHOEnergySimHits[detId] +=
simhits->energy();
605 fHFEnergySimHits[detId] +=
simhits->energy();
611 Double_t maxHBEnergy = 0.;
612 Double_t maxHEEnergy = 0.;
613 Double_t maxHFEnergy = 0.;
615 Double_t maxHBPhi = -1000.;
616 Double_t maxHEPhi = -1000.;
617 Double_t maxHOPhi = -1000.;
618 Double_t maxHFPhi = -1000.;
620 Double_t
PI = 3.141592653589;
625 std::vector<edm::Handle<HBHERecHitCollection>>
hbhe;
627 bool validHBHE =
true;
629 LogDebug(MsgLoggerCat) <<
"Unable to find any HBHERecHitCollections in event!";
634 std::vector<edm::Handle<HBHERecHitCollection>>::iterator ihbhe;
639 for (ihbhe =
hbhe.begin(); ihbhe !=
hbhe.end(); ++ihbhe) {
648 if ((jhbhe->energy()) > maxHBEnergy) {
649 maxHBEnergy = jhbhe->energy();
659 if ((jhbhe->energy()) > maxHEEnergy) {
660 maxHEEnergy = jhbhe->energy();
676 float deltaphi = maxHBPhi - fPhi;
677 if (fPhi > maxHBPhi) {
678 deltaphi = fPhi - maxHBPhi;
681 deltaphi = 2.0 *
PI - deltaphi;
694 float deltaphi = maxHEPhi - fPhi;
695 if (fPhi > maxHEPhi) {
696 deltaphi = fPhi - maxHEPhi;
699 deltaphi = 2.0 *
PI - deltaphi;
707 eventout +=
"\n Number of HBRecHits collected:............ ";
712 eventout +=
"\n Number of HERecHits collected:............ ";
722 std::vector<edm::Handle<HFRecHitCollection>>
hf;
726 LogDebug(MsgLoggerCat) <<
"Unable to find any HFRecHitCollections in event!";
730 std::vector<edm::Handle<HFRecHitCollection>>::iterator
ihf;
739 auto cellGeometry =
geometry->getSubdetectorGeometry(cell)->getGeometry(cell);
740 double fPhi = cellGeometry->getPosition().phi();
741 if ((jhf->energy()) > maxHFEnergy) {
742 maxHFEnergy = jhf->energy();
754 auto cellGeometry =
geometry->getSubdetectorGeometry(cell)->getGeometry(cell);
755 double fPhi = cellGeometry->getPosition().phi();
757 float deltaphi = maxHBPhi - fPhi;
758 if (fPhi > maxHFPhi) {
759 deltaphi = fPhi - maxHFPhi;
762 deltaphi = 2.0 *
PI - deltaphi;
771 eventout +=
"\n Number of HFDigis collected:.............. ";
780 std::vector<edm::Handle<HORecHitCollection>>
ho;
784 LogDebug(MsgLoggerCat) <<
"Unable to find any HORecHitCollections in event!";
789 std::vector<edm::Handle<HORecHitCollection>>::iterator iho;
792 for (iho =
ho.begin(); iho !=
ho.end(); ++iho) {
799 auto cellGeometry =
geometry->getSubdetectorGeometry(cell)->getGeometry(cell);
800 double fPhi = cellGeometry->getPosition().phi();
802 float deltaphi = maxHOPhi - fPhi;
803 if (fPhi > maxHOPhi) {
804 deltaphi = fPhi - maxHOPhi;
807 deltaphi = 2.0 *
PI - deltaphi;
815 eventout +=
"\n Number of HODigis collected:.............. ";
830 std::string MsgLoggerCat =
"GlobalRecHitsAnalyzer_fillTrk";
833 eventout =
"\nGathering info:";
838 bool validstrip =
true;
839 if (!rechitsmatched.
isValid()) {
840 LogDebug(MsgLoggerCat) <<
"Unable to find stripmatchedrechits in event!";
847 if (!tGeomHandle.isValid()) {
848 edm::LogWarning(MsgLoggerCat) <<
"Unable to find TrackerDigiGeometry in event!";
854 int nStripBrl = 0, nStripFwd = 0;
857 for (TrackerGeometry::DetContainer::const_iterator it = tGeomHandle->dets().begin();
858 it != tGeomHandle->dets().end();
860 uint32_t myid = ((*it)->geographicalId()).rawId();
861 DetId detid = ((*it)->geographicalId());
866 if (rechitmatchedMatch != rechitsmatched->
end()) {
869 rechitmatchedRange.
begin();
871 rechitmatchedRange.
end();
874 for (itermatched = rechitmatchedRangeIteratorBegin; itermatched != rechitmatchedRangeIteratorEnd;
879 float mindist = 999999.;
880 float distx = 999999.;
881 float disty = 999999.;
882 float dist = 999999.;
883 std::pair<LocalPoint, LocalVector> closestPair;
886 float rechitmatchedx =
position.x();
887 float rechitmatchedy =
position.y();
895 std::pair<LocalPoint, LocalVector> hitPair;
897 for (std::vector<PSimHit>::const_iterator
m =
matched.begin();
m !=
matched.end();
m++) {
900 distx = fabs(rechitmatchedx - hitPair.first.x());
901 disty = fabs(rechitmatchedy - hitPair.first.y());
902 dist =
sqrt(distx * distx + disty * disty);
904 if (dist < mindist) {
906 closestPair = hitPair;
1016 eventout +=
"\n Number of BrlStripRecHits collected:...... ";
1017 eventout += nStripBrl;
1020 for (
int i = 8;
i < 12; ++
i) {
1023 for (
int i = 16;
i < 19; ++
i) {
1028 eventout +=
"\n Number of FrwdStripRecHits collected:..... ";
1029 eventout += nStripFwd;
1031 for (
int i = 0;
i < 8; ++
i) {
1034 for (
int i = 12;
i < 16; ++
i) {
1043 bool validpixel =
true;
1045 LogDebug(MsgLoggerCat) <<
"Unable to find SiPixelRecHitCollection in event!";
1050 int nPxlBrl = 0, nPxlFwd = 0;
1052 for (TrackerGeometry::DetContainer::const_iterator it = tGeomHandle->dets().begin();
1053 it != tGeomHandle->dets().end();
1055 uint32_t myid = ((*it)->geographicalId()).rawId();
1056 DetId detId = ((*it)->geographicalId());
1063 if (pixeldet == recHitColl->
end())
1073 for (; pixeliter != pixelrechitRangeIteratorEnd; ++pixeliter) {
1078 float closest = 9999.9;
1080 float rechit_x = lp.
x();
1081 float rechit_y = lp.
y();
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);
1092 float sim_y1 = (*m).entryPoint().y();
1093 float sim_y2 = (*m).exitPoint().y();
1094 float sim_ypos = 0.5 * (sim_y1 + sim_y2);
1096 float x_res = fabs(sim_xpos - rechit_x);
1097 float y_res = fabs(sim_ypos - rechit_y);
1099 float dist =
sqrt(x_res * x_res + y_res * y_res);
1101 if (dist < closest) {
1130 if (tTopo->
pxfDisk(myid) == 1) {
1131 if (tTopo->
pxfSide(myid) == 1) {
1135 if (tTopo->
pxfSide(myid) == 2) {
1140 if (tTopo->
pxfDisk(myid) == 2) {
1141 if (tTopo->
pxfSide(myid) == 1) {
1145 if (tTopo->
pxfSide(myid) == 2) {
1156 eventout +=
"\n Number of BrlPixelRecHits collected:...... ";
1157 eventout += nPxlBrl;
1159 for (
int i = 0;
i < 3; ++
i) {
1164 eventout +=
"\n Number of FrwdPixelRecHits collected:..... ";
1165 eventout += nPxlFwd;
1168 for (
int i = 3;
i < 7; ++
i) {
1180 std::string MsgLoggerCat =
"GlobalRecHitsAnalyzer_fillMuon";
1184 eventout =
"\nGathering info:";
1188 if (!dtGeom.isValid()) {
1189 edm::LogWarning(MsgLoggerCat) <<
"Unable to find DTMuonGeometryRecord in event!";
1195 bool validdtsim =
true;
1197 LogDebug(MsgLoggerCat) <<
"Unable to find dtsimHits in event!";
1203 bool validdtrec =
true;
1205 LogDebug(MsgLoggerCat) <<
"Unable to find dtRecHits in event!";
1209 if (validdtsim && validdtrec) {
1210 std::map<DTWireId, edm::PSimHitContainer> simHitsPerWire =
1215 int nDt =
compute(dtGeom.product(), simHitsPerWire, recHitsPerWire, 1);
1218 eventout +=
"\n Number of DtMuonRecHits collected:........ ";
1230 bool validXFrame =
true;
1232 LogDebug(MsgLoggerCat) <<
"Unable to find muo CSC crossingFrame in event!";
1233 validXFrame =
false;
1239 for (
auto const& iHit :
simHits) {
1240 theMap[iHit.detUnitId()].push_back(iHit);
1246 if (!hGeom.isValid()) {
1247 edm::LogWarning(MsgLoggerCat) <<
"Unable to find CSCMuonGeometryRecord in event!";
1255 bool validCSC =
true;
1257 LogDebug(MsgLoggerCat) <<
"Unable to find CSC RecHits in event!";
1267 int detId = (*recHitItr).cscDetId().rawId();
1270 std::map<int, edm::PSimHitContainer>::const_iterator mapItr =
theMap.find(detId);
1271 if (mapItr !=
theMap.end()) {
1279 const CSCLayer*
layer = dynamic_cast<const CSCLayer*>(detUnit);
1281 int chamberType =
layer->chamber()->specs()->chamberType();
1287 eventout +=
"\n Number of CSCRecHits collected:........... ";
1294 std::map<double, int> mapsim, maprec;
1295 std::map<int, double> nmapsim, nmaprec;
1297 if (!rpcGeom.isValid()) {
1298 edm::LogWarning(MsgLoggerCat) <<
"Unable to find RPCMuonGeometryRecord in event!";
1304 bool validrpcsim =
true;
1306 LogDebug(MsgLoggerCat) <<
"Unable to find RPCSimHit in event!";
1307 validrpcsim =
false;
1312 bool validrpc =
true;
1314 LogDebug(MsgLoggerCat) <<
"Unable to find RPCRecHit in event!";
1322 for (recIt =
recHit->begin(); recIt !=
recHit->end(); ++recIt) {
1324 const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(Rid));
1327 eventout +=
"\n Number of RPCRecHits collected:........... ";
1336 LocalPoint rhitlocal = (*recIt).localPosition();
1337 double rhitlocalx = rhitlocal.
x();
1338 maprec[rhitlocalx] = nrec;
1342 for (std::map<double, int>::iterator iter = maprec.begin(); iter != maprec.end(); ++iter) {
1344 nmaprec[
i] = (*iter).first;
1349 edm::PSimHitContainer::const_iterator simIt;
1350 for (simIt =
simHit->begin(); simIt !=
simHit->end(); simIt++) {
1351 int ptype = (*simIt).particleType();
1354 LocalPoint shitlocal = (*simIt).localPosition();
1355 double shitlocalx = shitlocal.
x();
1356 mapsim[shitlocalx] = nsim;
1361 for (std::map<double, int>::iterator iter = mapsim.begin(); iter != mapsim.end(); ++iter) {
1363 nmapsim[
i] = (*iter).first;
1368 for (
int r = 0;
r < nsim;
r++) {
1375 eventout +=
"\n Number of RPCRecHits collected:........... ";
1393 LocalPoint localHit = plane.toLocal(globalpos);
1411 return std::pair<LocalPoint, LocalVector>(projectedPos, localStripDir);
1417 std::map<DTWireId, std::vector<DTRecHit1DPair>>
ret;
1421 ret[(*rechit).wireId()].push_back(*rechit);
1429 float xwire =
layer->specificTopology().wirePosition(wireId.
wire());
1432 float xEntry = entryP.
x() - xwire;
1433 float xExit = exitP.
x() - xwire;
1436 return fabs(xEntry - (entryP.
z() * (xExit - xEntry)) / (exitP.
z() - entryP.
z()));
1440 template <
typename type>
1443 const std::vector<type>&
recHits,
1444 const float simHitDist) {
1446 const type* theBestRecHit =
nullptr;
1450 if (fabs(distTmp - simHitDist) <
res) {
1451 res = fabs(distTmp - simHitDist);
1452 theBestRecHit = &(*recHit);
1456 return theBestRecHit;
1467 return fabs(
recHit.localPosition().x() -
layer->specificTopology().wirePosition(
recHit.wireId().wire()));
1470 template <
typename type>
1475 std::map<DTWireId, std::vector<PSimHit>> simHitsPerWire = _simHitsPerWire;
1476 std::map<DTWireId, std::vector<type>> recHitsPerWire = _recHitsPerWire;
1479 for (
std::map<
DTWireId, std::vector<PSimHit>>::const_iterator wireAndSHits = simHitsPerWire.begin();
1480 wireAndSHits != simHitsPerWire.end();
1482 DTWireId wireId = (*wireAndSHits).first;
1483 std::vector<PSimHit> simHitsInCell = (*wireAndSHits).second;
1490 if (muSimHit ==
nullptr) {
1497 if (simHitWireDist > 2.1) {
1502 if (recHitsPerWire.find(wireId) == recHitsPerWire.end()) {
1505 std::vector<type>
recHits = recHitsPerWire[wireId];