Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-05-15 22:20:56

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::Event &event,
0008                      const edm::EventSetup &setup,
0009                      const edm::ParameterSet &conf,
0010                      edm::ConsumesCollector &iC)
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   initEvent(event, setup);
0024 }
0025 
0026 MuonTruth::MuonTruth(const edm::ParameterSet &conf, edm::ConsumesCollector &&iC)
0027     : theDigiSimLinks(nullptr),
0028       theWireDigiSimLinks(nullptr),
0029       linksTag(conf.getParameter<edm::InputTag>("CSClinksTag")),
0030       wireLinksTag(conf.getParameter<edm::InputTag>("CSCwireLinksTag")),
0031       // CrossingFrame used or not ?
0032       crossingframe(conf.getParameter<bool>("crossingframe")),
0033       CSCsimHitsTag(conf.getParameter<edm::InputTag>("CSCsimHitsTag")),
0034       CSCsimHitsXFTag(conf.getParameter<edm::InputTag>("CSCsimHitsXFTag")),
0035       linksToken_(iC.consumes<DigiSimLinks>(linksTag)),
0036       wireLinksToken_(iC.consumes<DigiSimLinks>(wireLinksTag)) {
0037   if (crossingframe) {
0038     simHitsXFToken_ = iC.consumes<CrossingFrame<PSimHit>>(CSCsimHitsXFTag);
0039   } else if (!CSCsimHitsTag.label().empty()) {
0040     simHitsToken_ = iC.consumes<edm::PSimHitContainer>(CSCsimHitsTag);
0041   }
0042 }
0043 
0044 void MuonTruth::initEvent(const edm::Event &event, const edm::EventSetup &setup) {
0045   LogTrace("MuonTruth") << "getting CSC Strip DigiSimLink collection - " << linksTag;
0046   const edm::Handle<DigiSimLinks> &digiSimLinks = event.getHandle(linksToken_);
0047   theDigiSimLinks = digiSimLinks.product();
0048 
0049   LogTrace("MuonTruth") << "getting CSC Wire DigiSimLink collection - " << wireLinksTag;
0050   const edm::Handle<DigiSimLinks> &wireDigiSimLinks = event.getHandle(wireLinksToken_);
0051   theWireDigiSimLinks = wireDigiSimLinks.product();
0052 
0053   // get CSC Geometry to use CSCLayer methods
0054   const edm::ESHandle<CSCGeometry> &mugeom = setup.getHandle(geomToken_);
0055   cscgeom = mugeom.product();
0056 
0057   // get CSC Bad Chambers (ME4/2)
0058   const edm::ESHandle<CSCBadChambers> &badChambers = setup.getHandle(badToken_);
0059   cscBadChambers = badChambers.product();
0060 
0061   theSimHitMap.clear();
0062 
0063   if (crossingframe) {
0064     LogTrace("MuonTruth") << "getting CrossingFrame<PSimHit> collection - " << CSCsimHitsXFTag;
0065     const edm::Handle<CrossingFrame<PSimHit>> &cf = event.getHandle(simHitsXFToken_);
0066 
0067     std::unique_ptr<MixCollection<PSimHit>> CSCsimhits(new MixCollection<PSimHit>(cf.product()));
0068     LogTrace("MuonTruth") << "... size = " << CSCsimhits->size();
0069 
0070     for (MixCollection<PSimHit>::MixItr hitItr = CSCsimhits->begin(); hitItr != CSCsimhits->end(); ++hitItr) {
0071       theSimHitMap[hitItr->detUnitId()].push_back(*hitItr);
0072     }
0073 
0074   } else if (!CSCsimHitsTag.label().empty()) {
0075     LogTrace("MuonTruth") << "getting PSimHit collection - " << CSCsimHitsTag;
0076     const edm::Handle<edm::PSimHitContainer> &CSCsimhits = event.getHandle(simHitsToken_);
0077     LogTrace("MuonTruth") << "... size = " << CSCsimhits->size();
0078 
0079     for (edm::PSimHitContainer::const_iterator hitItr = CSCsimhits->begin(); hitItr != CSCsimhits->end(); ++hitItr) {
0080       theSimHitMap[hitItr->detUnitId()].push_back(*hitItr);
0081     }
0082   }
0083 }
0084 
0085 float MuonTruth::muonFraction() {
0086   if (theChargeMap.empty())
0087     return 0.;
0088 
0089   float muonCharge = 0.;
0090   for (std::map<SimHitIdpr, float>::const_iterator chargeMapItr = theChargeMap.begin();
0091        chargeMapItr != theChargeMap.end();
0092        ++chargeMapItr) {
0093     if (abs(particleType(chargeMapItr->first)) == 13) {
0094       muonCharge += chargeMapItr->second;
0095     }
0096   }
0097 
0098   return muonCharge / theTotalCharge;
0099 }
0100 
0101 std::vector<PSimHit> MuonTruth::simHits() {
0102   std::vector<PSimHit> result;
0103   for (std::map<SimHitIdpr, float>::const_iterator chargeMapItr = theChargeMap.begin();
0104        chargeMapItr != theChargeMap.end();
0105        ++chargeMapItr) {
0106     std::vector<PSimHit> trackHits = hitsFromSimTrack(chargeMapItr->first);
0107     result.insert(result.end(), trackHits.begin(), trackHits.end());
0108   }
0109 
0110   return result;
0111 }
0112 
0113 std::vector<PSimHit> MuonTruth::muonHits() {
0114   std::vector<PSimHit> result;
0115   std::vector<PSimHit> allHits = simHits();
0116   std::vector<PSimHit>::const_iterator hitItr = allHits.begin(), lastHit = allHits.end();
0117 
0118   for (; hitItr != lastHit; ++hitItr) {
0119     if (abs((*hitItr).particleType()) == 13) {
0120       result.push_back(*hitItr);
0121     }
0122   }
0123   return result;
0124 }
0125 
0126 std::vector<PSimHit> MuonTruth::hitsFromSimTrack(MuonTruth::SimHitIdpr truthId) {
0127   std::vector<PSimHit> result;
0128 
0129   auto found = theSimHitMap.find(theDetId);
0130   if (found != theSimHitMap.end()) {
0131     for (auto const &hit : found->second) {
0132       unsigned int hitTrack = hit.trackId();
0133       EncodedEventId hitEvId = hit.eventId();
0134 
0135       if (hitTrack == truthId.first && hitEvId == truthId.second) {
0136         result.push_back(hit);
0137       }
0138     }
0139   }
0140   return result;
0141 }
0142 
0143 int MuonTruth::particleType(MuonTruth::SimHitIdpr truthId) {
0144   int result = 0;
0145   const std::vector<PSimHit> &hits = hitsFromSimTrack(truthId);
0146   if (!hits.empty()) {
0147     result = hits[0].particleType();
0148   }
0149   return result;
0150 }
0151 
0152 void MuonTruth::analyze(const CSCRecHit2D &recHit) {
0153   theChargeMap.clear();
0154   theTotalCharge = 0.;
0155   theDetId = recHit.geographicalId().rawId();
0156 
0157   int nchannels = recHit.nStrips();
0158   const CSCLayerGeometry *laygeom = cscgeom->layer(recHit.cscDetId())->geometry();
0159 
0160   for (int idigi = 0; idigi < nchannels; ++idigi) {
0161     // strip and readout channel numbers may differ in ME1/1A
0162     int istrip = recHit.channels(idigi);
0163     int channel = laygeom->channel(istrip);
0164     float weight = recHit.adcs(idigi, 0);  // DL: I think this is wrong before and after...seems to
0165                                            // assume one time binadcContainer[idigi];
0166 
0167     DigiSimLinks::const_iterator layerLinks = theDigiSimLinks->find(theDetId);
0168 
0169     if (layerLinks != theDigiSimLinks->end()) {
0170       addChannel(*layerLinks, channel, weight);
0171     }
0172   }
0173 }
0174 
0175 void MuonTruth::analyze(const CSCStripDigi &stripDigi, int rawDetIdCorrespondingToCSCLayer) {
0176   theDetId = rawDetIdCorrespondingToCSCLayer;
0177   theChargeMap.clear();
0178   theTotalCharge = 0.;
0179 
0180   DigiSimLinks::const_iterator layerLinks = theDigiSimLinks->find(theDetId);
0181   if (layerLinks != theDigiSimLinks->end()) {
0182     addChannel(*layerLinks, stripDigi.getStrip(), 1.);
0183   }
0184 }
0185 
0186 void MuonTruth::analyze(const CSCWireDigi &wireDigi, int rawDetIdCorrespondingToCSCLayer) {
0187   theDetId = rawDetIdCorrespondingToCSCLayer;
0188   theChargeMap.clear();
0189   theTotalCharge = 0.;
0190 
0191   WireDigiSimLinks::const_iterator layerLinks = theWireDigiSimLinks->find(theDetId);
0192 
0193   if (layerLinks != theDigiSimLinks->end()) {
0194     // In the simulation digis, the channel labels for wires and strips must be
0195     // distinct, therefore:
0196     int wireDigiInSimulation = wireDigi.getWireGroup() + 100;
0197     //
0198     addChannel(*layerLinks, wireDigiInSimulation, 1.);
0199   }
0200 }
0201 
0202 void MuonTruth::addChannel(const LayerLinks &layerLinks, int channel, float weight) {
0203   LayerLinks::const_iterator linkItr = layerLinks.begin(), lastLayerLink = layerLinks.end();
0204 
0205   for (; linkItr != lastLayerLink; ++linkItr) {
0206     int linkChannel = linkItr->channel();
0207     if (linkChannel == channel) {
0208       float charge = linkItr->fraction() * weight;
0209       theTotalCharge += charge;
0210       // see if it's in the map
0211       SimHitIdpr truthId(linkItr->SimTrackId(), linkItr->eventId());
0212       std::map<SimHitIdpr, float>::const_iterator chargeMapItr = theChargeMap.find(truthId);
0213       if (chargeMapItr == theChargeMap.end()) {
0214         theChargeMap[truthId] = charge;
0215       } else {
0216         theChargeMap[truthId] += charge;
0217       }
0218     }
0219   }
0220 }