21 getAllProvenances(
false),
22 printProvenanceInfo(
false),
23 trackerHitAssociatorConfig_(iPSet, consumesCollector()),
25 consumesMany<edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>>();
26 consumesMany<edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>>();
27 consumesMany<edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>>();
28 std::string MsgLoggerCat =
"GlobalRecHitsAnalyzer_GlobalRecHitsAnalyzer";
90 edm::LogInfo(MsgLoggerCat) <<
"\n===============================\n"
91 <<
"Initialized as EDProducer with parameter values:\n"
92 <<
" Name = " <<
fName <<
"\n"
116 <<
"===============================\n";
124 string SiStripString[19] = {
"TECW1",
143 for (
int i = 0;
i < 19; ++
i) {
148 string hcharname, hchartitle;
150 for (
int amend = 0; amend < 19; ++amend) {
151 hcharname =
"hSiStripn_" + SiStripString[amend];
152 hchartitle = SiStripString[amend] +
" rechits";
156 hcharname =
"hSiStripResX_" + SiStripString[amend];
157 hchartitle = SiStripString[amend] +
" rechit x resolution";
161 hcharname =
"hSiStripResY_" + SiStripString[amend];
162 hchartitle = SiStripString[amend] +
" rechit y resolution";
170 string HCalString[4] = {
"HB",
"HE",
"HF",
"HO"};
171 float HCalnUpper[4] = {3000., 3000., 3000., 3000.};
172 float HCalnLower[4] = {0., 0., 0., 0.};
173 for (
int j = 0;
j < 4; ++
j) {
179 for (
int amend = 0; amend < 4; ++amend) {
180 hcharname =
"hHcaln_" + HCalString[amend];
181 hchartitle = HCalString[amend] +
" rechits";
182 mehHcaln[amend] = iBooker.
book1D(hcharname, hchartitle, 1000, HCalnLower[amend], HCalnUpper[amend]);
185 hcharname =
"hHcalRes_" + HCalString[amend];
186 hchartitle = HCalString[amend] +
" rechit resolution";
193 string ECalString[3] = {
"EB",
"EE",
"ES"};
194 int ECalnBins[3] = {1000, 3000, 150};
195 float ECalnUpper[3] = {20000., 62000., 3000.};
196 float ECalnLower[3] = {0., 0., 0.};
197 int ECalResBins[3] = {200, 200, 200};
198 float ECalResUpper[3] = {1., 0.3, .0002};
199 float ECalResLower[3] = {-1., -0.3, -.0002};
200 for (
int i = 0;
i < 3; ++
i) {
206 for (
int amend = 0; amend < 3; ++amend) {
207 hcharname =
"hEcaln_" + ECalString[amend];
208 hchartitle = ECalString[amend] +
" rechits";
209 mehEcaln[amend] = iBooker.
book1D(hcharname, hchartitle, ECalnBins[amend], ECalnLower[amend], ECalnUpper[amend]);
212 hcharname =
"hEcalRes_" + ECalString[amend];
213 hchartitle = ECalString[amend] +
" rechit resolution";
215 iBooker.
book1D(hcharname, hchartitle, ECalResBins[amend], ECalResLower[amend], ECalResUpper[amend]);
221 string SiPixelString[7] = {
"BRL1",
"BRL2",
"BRL3",
"FWD1n",
"FWD1p",
"FWD2n",
"FWD2p"};
222 for (
int j = 0;
j < 7; ++
j) {
229 for (
int amend = 0; amend < 7; ++amend) {
230 hcharname =
"hSiPixeln_" + SiPixelString[amend];
231 hchartitle = SiPixelString[amend] +
" rechits";
235 hcharname =
"hSiPixelResX_" + SiPixelString[amend];
236 hchartitle = SiPixelString[amend] +
" rechit x resolution";
240 hcharname =
"hSiPixelResY_" + SiPixelString[amend];
241 hchartitle = SiPixelString[amend] +
" rechit y resolution";
255 string n_List[3] = {
"hDtMuonn",
"hCSCn",
"hRPCn"};
256 string hist_string[3] = {
"Dt",
"CSC",
"RPC"};
258 for (
int amend = 0; amend < 3; ++amend) {
259 hchartitle = hist_string[amend] +
" rechits";
266 mehCSCn = iBooker.
book1D(n_List[amend], hchartitle, 50, 0., 500.);
271 mehRPCn = iBooker.
book1D(n_List[amend], hchartitle, 50, 0., 500.);
281 hcharname =
"hDtMuonRes";
282 hchartitle =
"DT wire distance resolution";
284 hcharname =
"CSCResRDPhi";
285 hchartitle =
"CSC perp*dphi resolution";
287 hcharname =
"hRPCResX";
288 hchartitle =
"RPC rechits x resolution";
293 std::string MsgLoggerCat =
"GlobalRecHitsAnalyzer_analyze";
303 edm::LogInfo(MsgLoggerCat) <<
"Processing run " << nrun <<
", event " <<
nevt <<
" (" <<
count <<
" events total)";
313 std::vector<const edm::StableProvenance*> AllProv;
314 iEvent.getAllStableProvenance(AllProv);
317 edm::LogInfo(MsgLoggerCat) <<
"Number of Provenances = " << AllProv.size();
320 TString eventout(
"\nProvenance info:\n");
322 for (
unsigned int i = 0;
i < AllProv.size(); ++
i) {
323 eventout +=
"\n ******************************";
324 eventout +=
"\n Module : ";
325 eventout += AllProv[
i]->moduleLabel();
326 eventout +=
"\n ProductID : ";
327 eventout += AllProv[
i]->productID().id();
328 eventout +=
"\n ClassName : ";
329 eventout += AllProv[
i]->className();
330 eventout +=
"\n InstanceName : ";
331 eventout += AllProv[
i]->productInstanceName();
332 eventout +=
"\n BranchName : ";
333 eventout += AllProv[
i]->branchName();
335 eventout +=
"\n ******************************\n";
353 edm::LogInfo(MsgLoggerCat) <<
"Done gathering data from event.";
359 std::string MsgLoggerCat =
"GlobalRecHitsAnalyzer_fillECal";
363 eventout =
"\nGathering info:";
373 bool validUncalibRecHitEB =
true;
374 if (!EcalUncalibRecHitEB.
isValid()) {
375 LogDebug(MsgLoggerCat) <<
"Unable to find EcalUncalRecHitEB in event!";
376 validUncalibRecHitEB =
false;
381 bool validRecHitEB =
true;
383 LogDebug(MsgLoggerCat) <<
"Unable to find EcalRecHitEB in event!";
384 validRecHitEB =
false;
389 bool validXFrame =
true;
390 if (!crossingFrame.
isValid()) {
391 LogDebug(MsgLoggerCat) <<
"Unable to find cal barrel crossingFrame in event!";
402 uint32_t crystid = ebid.
rawId();
403 ebSimMap[crystid] += iHit.energy();
409 if (validUncalibRecHitEB && validRecHitEB) {
427 eventout +=
"\n Number of EBRecHits collected:............ ";
428 eventout += nEBRecHits;
438 bool validuncalibRecHitEE =
true;
439 if (!EcalUncalibRecHitEE.
isValid()) {
440 LogDebug(MsgLoggerCat) <<
"Unable to find EcalUncalRecHitEE in event!";
441 validuncalibRecHitEE =
false;
446 bool validRecHitEE =
true;
448 LogDebug(MsgLoggerCat) <<
"Unable to find EcalRecHitEE in event!";
449 validRecHitEE =
false;
455 if (!crossingFrame.
isValid()) {
456 LogDebug(MsgLoggerCat) <<
"Unable to find cal endcap crossingFrame in event!";
467 uint32_t crystid = eeid.
rawId();
468 eeSimMap[crystid] += iHit.energy();
473 if (validuncalibRecHitEE && validRecHitEE) {
492 eventout +=
"\n Number of EERecHits collected:............ ";
493 eventout += nEERecHits;
503 bool validRecHitES =
true;
505 LogDebug(MsgLoggerCat) <<
"Unable to find EcalRecHitES in event!";
506 validRecHitES =
false;
512 if (!crossingFrame.
isValid()) {
513 LogDebug(MsgLoggerCat) <<
"Unable to find cal preshower crossingFrame in event!";
524 uint32_t crystid = esid.
rawId();
525 esSimMap[crystid] += iHit.energy();
541 eventout +=
"\n Number of ESRecHits collected:............ ";
542 eventout += nESRecHits;
554 std::string MsgLoggerCat =
"GlobalRecHitsAnalyzer_fillHCal";
558 eventout =
"\nGathering info:";
564 edm::LogWarning(MsgLoggerCat) <<
"Unable to find CaloGeometry in event!";
573 bool validhcalHits =
true;
575 LogDebug(MsgLoggerCat) <<
"Unable to find hcalHits in event!";
576 validhcalHits =
false;
579 std::map<HcalDetId, float> fHBEnergySimHits;
580 std::map<HcalDetId, float> fHEEnergySimHits;
581 std::map<HcalDetId, float> fHOEnergySimHits;
582 std::map<HcalDetId, float> fHFEnergySimHits;
586 for (std::vector<PCaloHit>::const_iterator
simhits = simhitResult->begin();
simhits != simhitResult->end();
591 fHBEnergySimHits[detId] +=
simhits->energy();
594 fHEEnergySimHits[detId] +=
simhits->energy();
597 fHOEnergySimHits[detId] +=
simhits->energy();
600 fHFEnergySimHits[detId] +=
simhits->energy();
606 Double_t maxHBEnergy = 0.;
607 Double_t maxHEEnergy = 0.;
608 Double_t maxHFEnergy = 0.;
610 Double_t maxHBPhi = -1000.;
611 Double_t maxHEPhi = -1000.;
612 Double_t maxHOPhi = -1000.;
613 Double_t maxHFPhi = -1000.;
615 Double_t
PI = 3.141592653589;
620 std::vector<edm::Handle<HBHERecHitCollection>>
hbhe;
622 bool validHBHE =
true;
624 LogDebug(MsgLoggerCat) <<
"Unable to find any HBHERecHitCollections in event!";
629 std::vector<edm::Handle<HBHERecHitCollection>>::iterator ihbhe;
634 for (ihbhe =
hbhe.begin(); ihbhe !=
hbhe.end(); ++ihbhe) {
643 if ((jhbhe->energy()) > maxHBEnergy) {
644 maxHBEnergy = jhbhe->energy();
654 if ((jhbhe->energy()) > maxHEEnergy) {
655 maxHEEnergy = jhbhe->energy();
671 float deltaphi = maxHBPhi - fPhi;
672 if (fPhi > maxHBPhi) {
673 deltaphi = fPhi - maxHBPhi;
676 deltaphi = 2.0 *
PI - deltaphi;
689 float deltaphi = maxHEPhi - fPhi;
690 if (fPhi > maxHEPhi) {
691 deltaphi = fPhi - maxHEPhi;
694 deltaphi = 2.0 *
PI - deltaphi;
702 eventout +=
"\n Number of HBRecHits collected:............ ";
707 eventout +=
"\n Number of HERecHits collected:............ ";
717 std::vector<edm::Handle<HFRecHitCollection>>
hf;
721 LogDebug(MsgLoggerCat) <<
"Unable to find any HFRecHitCollections in event!";
725 std::vector<edm::Handle<HFRecHitCollection>>::iterator
ihf;
734 auto cellGeometry =
geometry->getSubdetectorGeometry(cell)->getGeometry(cell);
735 double fPhi = cellGeometry->getPosition().phi();
736 if ((jhf->energy()) > maxHFEnergy) {
737 maxHFEnergy = jhf->energy();
749 auto cellGeometry =
geometry->getSubdetectorGeometry(cell)->getGeometry(cell);
750 double fPhi = cellGeometry->getPosition().phi();
752 float deltaphi = maxHBPhi - fPhi;
753 if (fPhi > maxHFPhi) {
754 deltaphi = fPhi - maxHFPhi;
757 deltaphi = 2.0 *
PI - deltaphi;
766 eventout +=
"\n Number of HFDigis collected:.............. ";
775 std::vector<edm::Handle<HORecHitCollection>>
ho;
779 LogDebug(MsgLoggerCat) <<
"Unable to find any HORecHitCollections in event!";
784 std::vector<edm::Handle<HORecHitCollection>>::iterator iho;
787 for (iho =
ho.begin(); iho !=
ho.end(); ++iho) {
794 auto cellGeometry =
geometry->getSubdetectorGeometry(cell)->getGeometry(cell);
795 double fPhi = cellGeometry->getPosition().phi();
797 float deltaphi = maxHOPhi - fPhi;
798 if (fPhi > maxHOPhi) {
799 deltaphi = fPhi - maxHOPhi;
802 deltaphi = 2.0 *
PI - deltaphi;
810 eventout +=
"\n Number of HODigis collected:.............. ";
828 std::string MsgLoggerCat =
"GlobalRecHitsAnalyzer_fillTrk";
832 eventout =
"\nGathering info:";
837 bool validstrip =
true;
838 if (!rechitsmatched.
isValid()) {
839 LogDebug(MsgLoggerCat) <<
"Unable to find stripmatchedrechits in event!";
848 edm::LogWarning(MsgLoggerCat) <<
"Unable to find TrackerDigiGeometry in event!";
854 int nStripBrl = 0, nStripFwd = 0;
857 for (TrackerGeometry::DetContainer::const_iterator it = pDD->
dets().begin(); it != pDD->
dets().end(); ++it) {
858 uint32_t myid = ((*it)->geographicalId()).rawId();
859 DetId detid = ((*it)->geographicalId());
864 if (rechitmatchedMatch != rechitsmatched->
end()) {
867 rechitmatchedRange.
begin();
869 rechitmatchedRange.
end();
872 for (itermatched = rechitmatchedRangeIteratorBegin; itermatched != rechitmatchedRangeIteratorEnd;
877 float mindist = 999999.;
878 float distx = 999999.;
879 float disty = 999999.;
880 float dist = 999999.;
881 std::pair<LocalPoint, LocalVector> closestPair;
884 float rechitmatchedx =
position.x();
885 float rechitmatchedy =
position.y();
893 std::pair<LocalPoint, LocalVector> hitPair;
895 for (std::vector<PSimHit>::const_iterator
m =
matched.begin();
m !=
matched.end();
m++) {
898 distx = fabs(rechitmatchedx - hitPair.first.x());
899 disty = fabs(rechitmatchedy - hitPair.first.y());
900 dist =
sqrt(distx * distx + disty * disty);
902 if (dist < mindist) {
904 closestPair = hitPair;
1014 eventout +=
"\n Number of BrlStripRecHits collected:...... ";
1015 eventout += nStripBrl;
1018 for (
int i = 8;
i < 12; ++
i) {
1021 for (
int i = 16;
i < 19; ++
i) {
1026 eventout +=
"\n Number of FrwdStripRecHits collected:..... ";
1027 eventout += nStripFwd;
1029 for (
int i = 0;
i < 8; ++
i) {
1032 for (
int i = 12;
i < 16; ++
i) {
1041 bool validpixel =
true;
1043 LogDebug(MsgLoggerCat) <<
"Unable to find SiPixelRecHitCollection in event!";
1050 if (!
geom.isValid()) {
1051 edm::LogWarning(MsgLoggerCat) <<
"Unable to find TrackerDigiGeometry in event!";
1056 int nPxlBrl = 0, nPxlFwd = 0;
1058 for (TrackerGeometry::DetContainer::const_iterator it =
geom->dets().begin(); it !=
geom->dets().end(); ++it) {
1059 uint32_t myid = ((*it)->geographicalId()).rawId();
1060 DetId detId = ((*it)->geographicalId());
1067 if (pixeldet == recHitColl->
end())
1077 for (; pixeliter != pixelrechitRangeIteratorEnd; ++pixeliter) {
1082 float closest = 9999.9;
1084 float rechit_x = lp.
x();
1085 float rechit_y = lp.
y();
1091 for (std::vector<PSimHit>::const_iterator
m =
matched.begin();
m !=
matched.end(); ++
m) {
1092 float sim_x1 = (*m).entryPoint().x();
1093 float sim_x2 = (*m).exitPoint().x();
1094 float sim_xpos = 0.5 * (sim_x1 + sim_x2);
1096 float sim_y1 = (*m).entryPoint().y();
1097 float sim_y2 = (*m).exitPoint().y();
1098 float sim_ypos = 0.5 * (sim_y1 + sim_y2);
1100 float x_res = fabs(sim_xpos - rechit_x);
1101 float y_res = fabs(sim_ypos - rechit_y);
1103 float dist =
sqrt(x_res * x_res + y_res * y_res);
1105 if (dist < closest) {
1134 if (tTopo->
pxfDisk(myid) == 1) {
1135 if (tTopo->
pxfSide(myid) == 1) {
1139 if (tTopo->
pxfSide(myid) == 2) {
1144 if (tTopo->
pxfDisk(myid) == 2) {
1145 if (tTopo->
pxfSide(myid) == 1) {
1149 if (tTopo->
pxfSide(myid) == 2) {
1160 eventout +=
"\n Number of BrlPixelRecHits collected:...... ";
1161 eventout += nPxlBrl;
1163 for (
int i = 0;
i < 3; ++
i) {
1168 eventout +=
"\n Number of FrwdPixelRecHits collected:..... ";
1169 eventout += nPxlFwd;
1172 for (
int i = 3;
i < 7; ++
i) {
1184 std::string MsgLoggerCat =
"GlobalRecHitsAnalyzer_fillMuon";
1188 eventout =
"\nGathering info:";
1194 edm::LogWarning(MsgLoggerCat) <<
"Unable to find DTMuonGeometryRecord in event!";
1200 bool validdtsim =
true;
1202 LogDebug(MsgLoggerCat) <<
"Unable to find dtsimHits in event!";
1208 bool validdtrec =
true;
1210 LogDebug(MsgLoggerCat) <<
"Unable to find dtRecHits in event!";
1214 if (validdtsim && validdtrec) {
1215 std::map<DTWireId, edm::PSimHitContainer> simHitsPerWire =
1220 int nDt =
compute(dtGeom.
product(), simHitsPerWire, recHitsPerWire, 1);
1223 eventout +=
"\n Number of DtMuonRecHits collected:........ ";
1235 bool validXFrame =
true;
1237 LogDebug(MsgLoggerCat) <<
"Unable to find muo CSC crossingFrame in event!";
1238 validXFrame =
false;
1244 for (
auto const& iHit :
simHits) {
1245 theMap[iHit.detUnitId()].push_back(iHit);
1253 edm::LogWarning(MsgLoggerCat) <<
"Unable to find CSCMuonGeometryRecord in event!";
1261 bool validCSC =
true;
1263 LogDebug(MsgLoggerCat) <<
"Unable to find CSC RecHits in event!";
1273 int detId = (*recHitItr).cscDetId().rawId();
1276 std::map<int, edm::PSimHitContainer>::const_iterator mapItr =
theMap.find(detId);
1277 if (mapItr !=
theMap.end()) {
1285 const CSCLayer*
layer = dynamic_cast<const CSCLayer*>(detUnit);
1287 int chamberType =
layer->chamber()->specs()->chamberType();
1293 eventout +=
"\n Number of CSCRecHits collected:........... ";
1300 std::map<double, int> mapsim, maprec;
1301 std::map<int, double> nmapsim, nmaprec;
1306 edm::LogWarning(MsgLoggerCat) <<
"Unable to find RPCMuonGeometryRecord in event!";
1312 bool validrpcsim =
true;
1314 LogDebug(MsgLoggerCat) <<
"Unable to find RPCSimHit in event!";
1315 validrpcsim =
false;
1320 bool validrpc =
true;
1322 LogDebug(MsgLoggerCat) <<
"Unable to find RPCRecHit in event!";
1330 for (recIt =
recHit->begin(); recIt !=
recHit->end(); ++recIt) {
1332 const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->
roll(Rid));
1335 eventout +=
"\n Number of RPCRecHits collected:........... ";
1344 LocalPoint rhitlocal = (*recIt).localPosition();
1345 double rhitlocalx = rhitlocal.
x();
1346 maprec[rhitlocalx] = nrec;
1350 for (std::map<double, int>::iterator iter = maprec.begin(); iter != maprec.end(); ++iter) {
1352 nmaprec[
i] = (*iter).first;
1357 edm::PSimHitContainer::const_iterator simIt;
1358 for (simIt =
simHit->begin(); simIt !=
simHit->end(); simIt++) {
1359 int ptype = (*simIt).particleType();
1362 LocalPoint shitlocal = (*simIt).localPosition();
1363 double shitlocalx = shitlocal.
x();
1364 mapsim[shitlocalx] = nsim;
1369 for (std::map<double, int>::iterator iter = mapsim.begin(); iter != mapsim.end(); ++iter) {
1371 nmapsim[
i] = (*iter).first;
1376 for (
int r = 0;
r < nsim;
r++) {
1383 eventout +=
"\n Number of RPCRecHits collected:........... ";
1401 LocalPoint localHit = plane.toLocal(globalpos);
1419 return std::pair<LocalPoint, LocalVector>(projectedPos, localStripDir);
1425 std::map<DTWireId, std::vector<DTRecHit1DPair>>
ret;
1429 ret[(*rechit).wireId()].push_back(*rechit);
1437 float xwire =
layer->specificTopology().wirePosition(wireId.
wire());
1440 float xEntry = entryP.
x() - xwire;
1441 float xExit = exitP.
x() - xwire;
1444 return fabs(xEntry - (entryP.
z() * (xExit - xEntry)) / (exitP.
z() - entryP.
z()));
1448 template <
typename type>
1451 const std::vector<type>&
recHits,
1452 const float simHitDist) {
1454 const type* theBestRecHit =
nullptr;
1458 if (fabs(distTmp - simHitDist) <
res) {
1459 res = fabs(distTmp - simHitDist);
1460 theBestRecHit = &(*recHit);
1464 return theBestRecHit;
1475 return fabs(
recHit.localPosition().x() -
layer->specificTopology().wirePosition(
recHit.wireId().wire()));
1478 template <
typename type>
1483 std::map<DTWireId, std::vector<PSimHit>> simHitsPerWire = _simHitsPerWire;
1484 std::map<DTWireId, std::vector<type>> recHitsPerWire = _recHitsPerWire;
1487 for (
std::map<
DTWireId, std::vector<PSimHit>>::const_iterator wireAndSHits = simHitsPerWire.begin();
1488 wireAndSHits != simHitsPerWire.end();
1490 DTWireId wireId = (*wireAndSHits).first;
1491 std::vector<PSimHit> simHitsInCell = (*wireAndSHits).second;
1498 if (muSimHit ==
nullptr) {
1505 if (simHitWireDist > 2.1) {
1510 if (recHitsPerWire.find(wireId) == recHitsPerWire.end()) {
1513 std::vector<type>
recHits = recHitsPerWire[wireId];