File indexing completed on 2023-08-10 23:18:30
0001 #include "DataFormats/Common/interface/DetSetVector.h"
0002 #include "FWCore/Framework/interface/ESHandle.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h"
0005 #include "SimMuon/MCTruth/interface/MuonTruth.h"
0006
0007 MuonTruth::MuonTruth(const edm::ParameterSet &conf, edm::ConsumesCollector &&iC)
0008 : cscBadChambers(nullptr),
0009 theTotalCharge(0.0),
0010 theDetId(0),
0011 theDigiSimLinks(nullptr),
0012 theWireDigiSimLinks(nullptr),
0013 linksTag(conf.getParameter<edm::InputTag>("CSClinksTag")),
0014 wireLinksTag(conf.getParameter<edm::InputTag>("CSCwireLinksTag")),
0015
0016 crossingframe(conf.getParameter<bool>("crossingframe")),
0017 CSCsimHitsTag(conf.getParameter<edm::InputTag>("CSCsimHitsTag")),
0018 CSCsimHitsXFTag(conf.getParameter<edm::InputTag>("CSCsimHitsXFTag")),
0019 geomToken_(iC.esConsumes<CSCGeometry, MuonGeometryRecord>()),
0020 badToken_(iC.esConsumes<CSCBadChambers, CSCBadChambersRcd>()),
0021 linksToken_(iC.consumes<DigiSimLinks>(linksTag)),
0022 wireLinksToken_(iC.consumes<DigiSimLinks>(wireLinksTag)),
0023 cscgeom(nullptr) {
0024 if (crossingframe) {
0025 simHitsXFToken_ = iC.consumes<CrossingFrame<PSimHit>>(CSCsimHitsXFTag);
0026 } else if (!CSCsimHitsTag.label().empty()) {
0027 simHitsToken_ = iC.consumes<edm::PSimHitContainer>(CSCsimHitsTag);
0028 }
0029 }
0030
0031 void MuonTruth::initEvent(const edm::Event &event, const edm::EventSetup &setup) {
0032 theChargeMap.clear();
0033 theTotalCharge = 0.0;
0034 theDetId = 0;
0035
0036 LogTrace("MuonTruth") << "getting CSC Strip DigiSimLink collection - " << linksTag;
0037 const edm::Handle<DigiSimLinks> &digiSimLinks = event.getHandle(linksToken_);
0038 theDigiSimLinks = digiSimLinks.product();
0039
0040 LogTrace("MuonTruth") << "getting CSC Wire DigiSimLink collection - " << wireLinksTag;
0041 const edm::Handle<DigiSimLinks> &wireDigiSimLinks = event.getHandle(wireLinksToken_);
0042 theWireDigiSimLinks = wireDigiSimLinks.product();
0043
0044
0045 const edm::ESHandle<CSCGeometry> &mugeom = setup.getHandle(geomToken_);
0046 cscgeom = mugeom.product();
0047
0048
0049 const edm::ESHandle<CSCBadChambers> &badChambers = setup.getHandle(badToken_);
0050 cscBadChambers = badChambers.product();
0051
0052 theSimHitMap.clear();
0053
0054 if (crossingframe) {
0055 LogTrace("MuonTruth") << "getting CrossingFrame<PSimHit> collection - " << CSCsimHitsXFTag;
0056 const edm::Handle<CrossingFrame<PSimHit>> &cf = event.getHandle(simHitsXFToken_);
0057
0058 std::unique_ptr<MixCollection<PSimHit>> CSCsimhits(new MixCollection<PSimHit>(cf.product()));
0059 LogTrace("MuonTruth") << "... size = " << CSCsimhits->size();
0060
0061 for (MixCollection<PSimHit>::MixItr hitItr = CSCsimhits->begin(); hitItr != CSCsimhits->end(); ++hitItr) {
0062 theSimHitMap[hitItr->detUnitId()].push_back(*hitItr);
0063 }
0064
0065 } else if (!CSCsimHitsTag.label().empty()) {
0066 LogTrace("MuonTruth") << "getting PSimHit collection - " << CSCsimHitsTag;
0067 const edm::Handle<edm::PSimHitContainer> &CSCsimhits = event.getHandle(simHitsToken_);
0068 LogTrace("MuonTruth") << "... size = " << CSCsimhits->size();
0069
0070 for (edm::PSimHitContainer::const_iterator hitItr = CSCsimhits->begin(); hitItr != CSCsimhits->end(); ++hitItr) {
0071 theSimHitMap[hitItr->detUnitId()].push_back(*hitItr);
0072 }
0073 }
0074 }
0075
0076 float MuonTruth::muonFraction() {
0077 if (theChargeMap.empty())
0078 return 0.;
0079
0080 float muonCharge = 0.;
0081 for (std::map<SimHitIdpr, float>::const_iterator chargeMapItr = theChargeMap.begin();
0082 chargeMapItr != theChargeMap.end();
0083 ++chargeMapItr) {
0084 if (abs(particleType(chargeMapItr->first)) == 13) {
0085 muonCharge += chargeMapItr->second;
0086 }
0087 }
0088
0089 return muonCharge / theTotalCharge;
0090 }
0091
0092 std::vector<PSimHit> MuonTruth::simHits() {
0093 std::vector<PSimHit> result;
0094 for (std::map<SimHitIdpr, float>::const_iterator chargeMapItr = theChargeMap.begin();
0095 chargeMapItr != theChargeMap.end();
0096 ++chargeMapItr) {
0097 std::vector<PSimHit> trackHits = hitsFromSimTrack(chargeMapItr->first);
0098 result.insert(result.end(), trackHits.begin(), trackHits.end());
0099 }
0100
0101 return result;
0102 }
0103
0104 std::vector<PSimHit> MuonTruth::muonHits() {
0105 std::vector<PSimHit> result;
0106 std::vector<PSimHit> allHits = simHits();
0107 std::vector<PSimHit>::const_iterator hitItr = allHits.begin(), lastHit = allHits.end();
0108
0109 for (; hitItr != lastHit; ++hitItr) {
0110 if (abs((*hitItr).particleType()) == 13) {
0111 result.push_back(*hitItr);
0112 }
0113 }
0114 return result;
0115 }
0116
0117 std::vector<PSimHit> MuonTruth::hitsFromSimTrack(MuonTruth::SimHitIdpr truthId) {
0118 std::vector<PSimHit> result;
0119
0120 auto found = theSimHitMap.find(theDetId);
0121 if (found != theSimHitMap.end()) {
0122 for (auto const &hit : found->second) {
0123 unsigned int hitTrack = hit.trackId();
0124 EncodedEventId hitEvId = hit.eventId();
0125
0126 if (hitTrack == truthId.first && hitEvId == truthId.second) {
0127 result.push_back(hit);
0128 }
0129 }
0130 }
0131 return result;
0132 }
0133
0134 int MuonTruth::particleType(MuonTruth::SimHitIdpr truthId) {
0135 int result = 0;
0136 const std::vector<PSimHit> &hits = hitsFromSimTrack(truthId);
0137 if (!hits.empty()) {
0138 result = hits[0].particleType();
0139 }
0140 return result;
0141 }
0142
0143 void MuonTruth::analyze(const CSCRecHit2D &recHit) {
0144 theChargeMap.clear();
0145 theTotalCharge = 0.;
0146 theDetId = recHit.geographicalId().rawId();
0147
0148 int nchannels = recHit.nStrips();
0149 const CSCLayerGeometry *laygeom = cscgeom->layer(recHit.cscDetId())->geometry();
0150
0151 for (int idigi = 0; idigi < nchannels; ++idigi) {
0152
0153 int istrip = recHit.channels(idigi);
0154 int channel = laygeom->channel(istrip);
0155 float weight = recHit.adcs(idigi, 0);
0156
0157
0158 DigiSimLinks::const_iterator layerLinks = theDigiSimLinks->find(theDetId);
0159
0160 if (layerLinks != theDigiSimLinks->end()) {
0161 addChannel(*layerLinks, channel, weight);
0162 }
0163 }
0164 }
0165
0166 void MuonTruth::analyze(const CSCStripDigi &stripDigi, int rawDetIdCorrespondingToCSCLayer) {
0167 theDetId = rawDetIdCorrespondingToCSCLayer;
0168 theChargeMap.clear();
0169 theTotalCharge = 0.;
0170
0171 DigiSimLinks::const_iterator layerLinks = theDigiSimLinks->find(theDetId);
0172 if (layerLinks != theDigiSimLinks->end()) {
0173 addChannel(*layerLinks, stripDigi.getStrip(), 1.);
0174 }
0175 }
0176
0177 void MuonTruth::analyze(const CSCWireDigi &wireDigi, int rawDetIdCorrespondingToCSCLayer) {
0178 theDetId = rawDetIdCorrespondingToCSCLayer;
0179 theChargeMap.clear();
0180 theTotalCharge = 0.;
0181
0182 WireDigiSimLinks::const_iterator layerLinks = theWireDigiSimLinks->find(theDetId);
0183
0184 if (layerLinks != theDigiSimLinks->end()) {
0185
0186
0187 int wireDigiInSimulation = wireDigi.getWireGroup() + 100;
0188
0189 addChannel(*layerLinks, wireDigiInSimulation, 1.);
0190 }
0191 }
0192
0193 void MuonTruth::addChannel(const LayerLinks &layerLinks, int channel, float weight) {
0194 LayerLinks::const_iterator linkItr = layerLinks.begin(), lastLayerLink = layerLinks.end();
0195
0196 for (; linkItr != lastLayerLink; ++linkItr) {
0197 int linkChannel = linkItr->channel();
0198 if (linkChannel == channel) {
0199 float charge = linkItr->fraction() * weight;
0200 theTotalCharge += charge;
0201
0202 SimHitIdpr truthId(linkItr->SimTrackId(), linkItr->eventId());
0203 std::map<SimHitIdpr, float>::const_iterator chargeMapItr = theChargeMap.find(truthId);
0204 if (chargeMapItr == theChargeMap.end()) {
0205 theChargeMap[truthId] = charge;
0206 } else {
0207 theChargeMap[truthId] += charge;
0208 }
0209 }
0210 }
0211 }