MCTruth/src/MuonTruth.cc

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       // CrossingFrame used or not ?
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   // get CSC Geometry to use CSCLayer methods
0045   const edm::ESHandle<CSCGeometry> &mugeom = setup.getHandle(geomToken_);
0046   cscgeom = mugeom.product();
0047
0048   // get CSC Bad Chambers (ME4/2)
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     // strip and readout channel numbers may differ in ME1/1A
0153     int istrip = recHit.channels(idigi);
0154     int channel = laygeom->channel(istrip);
0155     float weight = recHit.adcs(idigi, 0);  // DL: I think this is wrong before and after...seems to
0156                                            // assume one time binadcContainer[idigi];
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     // In the simulation digis, the channel labels for wires and strips must be
0186     // distinct, therefore:
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       // see if it's in the map
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 }