Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:30:47 

0001 #include "DataFormats/DetId/interface/DetId.h"
0002 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "Geometry/DTGeometry/interface/DTLayer.h"
0006 #include "Geometry/DTGeometry/interface/DTTopology.h"
0007 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0008 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0009 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0010 #include "SimMuon/MCTruth/interface/DTHitAssociator.h"
0011 #include <sstream>
0012 #include <string>
0013 
0014 using namespace std;
0015 
0016 // Constructor
0017 DTHitAssociator::Config::Config(const edm::ParameterSet &conf, edm::ConsumesCollector iC)
0018     : DTsimhitsTag(conf.getParameter<edm::InputTag>("DTsimhitsTag")),
0019       DTsimhitsXFTag(conf.getParameter<edm::InputTag>("DTsimhitsXFTag")),
0020       DTdigiTag(conf.getParameter<edm::InputTag>("DTdigiTag")),
0021       DTdigisimlinkTag(conf.getParameter<edm::InputTag>("DTdigisimlinkTag")),
0022       DTrechitTag(conf.getParameter<edm::InputTag>("DTrechitTag")),
0023       geomToken(iC.esConsumes()),
0024       // nice printout of DT hits
0025       dumpDT(conf.getParameter<bool>("dumpDT")),
0026       // CrossingFrame used or not ?
0027       crossingframe(conf.getParameter<bool>("crossingframe")),
0028       // Event contain the DTDigiSimLink collection ?
0029       links_exist(conf.getParameter<bool>("links_exist")),
0030       // associatorByWire links to a RecHit all the "valid" SimHits on the same
0031       // DT wire
0032       associatorByWire(conf.getParameter<bool>("associatorByWire")) {
0033   if (crossingframe) {
0034     DTsimhitsXFToken = iC.consumes<CrossingFrame<PSimHit>>(DTsimhitsXFTag);
0035   } else if (!DTsimhitsTag.label().empty()) {
0036     DTsimhitsToken = iC.consumes<edm::PSimHitContainer>(DTsimhitsTag);
0037   }
0038   DTdigiToken = iC.consumes<DTDigiCollection>(DTdigiTag);
0039   DTdigisimlinkToken = iC.consumes<DTDigiSimLinkCollection>(DTdigisimlinkTag);
0040 
0041   if (dumpDT) {
0042     DTrechitToken = iC.consumes<DTRecHitCollection>(DTrechitTag);
0043   }
0044 
0045   if (!links_exist && !associatorByWire) {
0046     LogTrace("DTHitAssociator") << "*** WARNING in DTHitAssociator::DTHitAssociator: associatorByWire "
0047                                    "reset to TRUE !"
0048                                 << "    \t (missing DTDigiSimLinkCollection ?)";
0049     associatorByWire = true;
0050   }
0051 }
0052 
0053 DTHitAssociator::DTHitAssociator(const edm::Event &iEvent,
0054                                  const edm::EventSetup &iSetup,
0055                                  const Config &conf,
0056                                  bool printRtS)
0057     :  // input collection labels
0058       config_(conf),
0059       printRtS(true)
0060 
0061 {
0062   initEvent(iEvent, iSetup);
0063 }
0064 
0065 void DTHitAssociator::initEvent(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0066   LogTrace("DTHitAssociator") << "DTHitAssociator constructor: dumpDT = " << config_.dumpDT
0067                               << ", crossingframe = " << config_.crossingframe
0068                               << ", links_exist = " << config_.links_exist
0069                               << ", associatorByWire = " << config_.associatorByWire;
0070 
0071   // need DT Geometry to discard hits for which the drift time parametrisation
0072   // is not applicable
0073   edm::ESHandle<DTGeometry> muonGeom = iSetup.getHandle(config_.geomToken);
0074 
0075   // Get the DT SimHits from the event and map PSimHit by DTWireId
0076   mapOfSimHit.clear();
0077   bool takeHit(true);
0078 
0079   if (config_.crossingframe) {
0080     LogTrace("DTHitAssociator") << "getting CrossingFrame<PSimHit> collection - " << config_.DTsimhitsXFTag;
0081     CrossingFrame<PSimHit> const &xFrame = iEvent.get(config_.DTsimhitsXFToken);
0082     unique_ptr<MixCollection<PSimHit>> DTsimhits(new MixCollection<PSimHit>(&xFrame));
0083     LogTrace("DTHitAssociator") << "... size = " << DTsimhits->size();
0084     MixCollection<PSimHit>::MixItr isimhit;
0085     for (isimhit = DTsimhits->begin(); isimhit != DTsimhits->end(); isimhit++) {
0086       DTWireId wireid((*isimhit).detUnitId());
0087       takeHit = SimHitOK(muonGeom, *isimhit);
0088       mapOfSimHit[wireid].push_back(make_pair(*isimhit, takeHit));
0089     }
0090   } else if (!config_.DTsimhitsTag.label().empty()) {
0091     LogTrace("DTHitAssociator") << "getting PSimHit collection - " << config_.DTsimhitsTag;
0092     edm::PSimHitContainer const &DTsimhits = iEvent.get(config_.DTsimhitsToken);
0093     LogTrace("DTHitAssociator") << "... size = " << DTsimhits.size();
0094     edm::PSimHitContainer::const_iterator isimhit;
0095     for (isimhit = DTsimhits.begin(); isimhit != DTsimhits.end(); isimhit++) {
0096       DTWireId wireid((*isimhit).detUnitId());
0097       takeHit = SimHitOK(muonGeom, *isimhit);
0098       mapOfSimHit[wireid].push_back(make_pair(*isimhit, takeHit));
0099     }
0100   }
0101 
0102   // Get the DT Digi collection from the event
0103   mapOfDigi.clear();
0104   LogTrace("DTHitAssociator") << "getting DTDigi collection - " << config_.DTdigiTag;
0105   edm::Handle<DTDigiCollection> digis = iEvent.getHandle(config_.DTdigiToken);
0106 
0107   if (digis.isValid()) {
0108     // Map DTDigi by DTWireId
0109     for (DTDigiCollection::DigiRangeIterator detUnit = digis->begin(); detUnit != digis->end(); ++detUnit) {
0110       const DTLayerId &layerid = (*detUnit).first;
0111       const DTDigiCollection::Range &range = (*detUnit).second;
0112 
0113       DTDigiCollection::const_iterator digi;
0114       for (digi = range.first; digi != range.second; ++digi) {
0115         DTWireId wireid(layerid, (*digi).wire());
0116         mapOfDigi[wireid].push_back(*digi);
0117       }
0118     }
0119   } else {
0120     LogTrace("DTHitAssociator") << "... NO DTDigi collection found !";
0121   }
0122 
0123   mapOfLinks.clear();
0124   if (config_.links_exist) {
0125     // Get the DT DigiSimLink collection from the event and map DTDigiSimLink by
0126     // DTWireId
0127     LogTrace("DTHitAssociator") << "getting DTDigiSimLink collection - " << config_.DTdigisimlinkTag;
0128     DTDigiSimLinkCollection const &digisimlinks = iEvent.get(config_.DTdigisimlinkToken);
0129 
0130     for (DTDigiSimLinkCollection::DigiRangeIterator detUnit = digisimlinks.begin(); detUnit != digisimlinks.end();
0131          ++detUnit) {
0132       const DTLayerId &layerid = (*detUnit).first;
0133       const DTDigiSimLinkCollection::Range &range = (*detUnit).second;
0134 
0135       DTDigiSimLinkCollection::const_iterator link;
0136       for (link = range.first; link != range.second; ++link) {
0137         DTWireId wireid(layerid, (*link).wire());
0138         mapOfLinks[wireid].push_back(*link);
0139       }
0140     }
0141   }
0142 
0143   if (config_.dumpDT && printRtS) {
0144     // Get the DT rechits from the event
0145     LogTrace("DTHitAssociator") << "getting DTRecHit1DPair collection - " << config_.DTrechitTag;
0146     DTRecHitCollection const &DTrechits = iEvent.get(config_.DTrechitToken);
0147     LogTrace("DTHitAssociator") << "... size = " << DTrechits.size();
0148 
0149     // map DTRecHit1DPair by DTWireId
0150     mapOfRecHit.clear();
0151     DTRecHitCollection::const_iterator rechit;
0152     for (rechit = DTrechits.begin(); rechit != DTrechits.end(); ++rechit) {
0153       DTWireId wireid = (*rechit).wireId();
0154       mapOfRecHit[wireid].push_back(*rechit);
0155     }
0156 
0157     if (mapOfSimHit.end() != mapOfSimHit.begin()) {
0158       LogTrace("DTHitAssociator") << "\n *** Dump DT PSimHit's ***";
0159 
0160       int jwire = 0;
0161       int ihit = 0;
0162 
0163       for (SimHitMap::const_iterator mapIT = mapOfSimHit.begin(); mapIT != mapOfSimHit.end(); ++mapIT, jwire++) {
0164         DTWireId wireid = (*mapIT).first;
0165         for (vector<PSimHit_withFlag>::const_iterator hitIT = mapOfSimHit[wireid].begin();
0166              hitIT != mapOfSimHit[wireid].end();
0167              hitIT++, ihit++) {
0168           PSimHit hit = hitIT->first;
0169           LogTrace("DTHitAssociator") << "PSimHit " << ihit << ", wire " << wireid  //<<", detID = "<<hit.detUnitId()
0170                                       << ", SimTrack Id:" << hit.trackId() << "/Evt:(" << hit.eventId().event() << ","
0171                                       << hit.eventId().bunchCrossing() << ") "
0172                                       << ", pdg = " << hit.particleType() << ", procs = " << hit.processType();
0173         }
0174       }
0175     } else {
0176       LogTrace("DTHitAssociator") << "\n *** There are NO DT PSimHit's ***";
0177     }
0178 
0179     if (mapOfDigi.end() != mapOfDigi.begin()) {
0180       int jwire = 0;
0181       int ihit = 0;
0182 
0183       for (DigiMap::const_iterator mapIT = mapOfDigi.begin(); mapIT != mapOfDigi.end(); ++mapIT, jwire++) {
0184         LogTrace("DTHitAssociator") << "\n *** Dump DT digis ***";
0185 
0186         DTWireId wireid = (*mapIT).first;
0187         for (vector<DTDigi>::const_iterator hitIT = mapOfDigi[wireid].begin(); hitIT != mapOfDigi[wireid].end();
0188              hitIT++, ihit++) {
0189           LogTrace("DTHitAssociator") << "DTDigi " << ihit << ", wire " << wireid << ", number = " << hitIT->number()
0190                                       << ", TDC counts = " << hitIT->countsTDC();
0191         }
0192       }
0193     } else {
0194       LogTrace("DTHitAssociator") << "\n *** There are NO DTDigi's ***";
0195     }
0196 
0197     if (mapOfRecHit.end() != mapOfRecHit.begin())
0198       LogTrace("DTHitAssociator") << "\n *** Analyze DTRecHitCollection by DTWireId ***";
0199 
0200     int iwire = 0;
0201     for (RecHitMap::const_iterator mapIT = mapOfRecHit.begin(); mapIT != mapOfRecHit.end(); ++mapIT, iwire++) {
0202       DTWireId wireid = (*mapIT).first;
0203       LogTrace("DTHitAssociator") << "\n======================================="
0204                                      "============================";
0205       LogTrace("DTHitAssociator") << "wire index = " << iwire << "  *** DTWireId = "
0206                                   << " (" << wireid << ")";
0207 
0208       if (mapOfSimHit.find(wireid) != mapOfSimHit.end()) {
0209         LogTrace("DTHitAssociator") << "\n" << mapOfSimHit[wireid].size() << " SimHits (PSimHit):";
0210 
0211         for (vector<PSimHit_withFlag>::const_iterator hitIT = mapOfSimHit[wireid].begin();
0212              hitIT != mapOfSimHit[wireid].end();
0213              ++hitIT) {
0214           stringstream tId;
0215           tId << (hitIT->first).trackId();
0216           string simhitlog = "\t SimTrack Id = " + tId.str();
0217           if (hitIT->second)
0218             simhitlog = simhitlog + "\t -VALID HIT-";
0219           else
0220             simhitlog = simhitlog + "\t -not valid hit-";
0221           LogTrace("DTHitAssociator") << simhitlog;
0222         }
0223       }
0224 
0225       if (mapOfLinks.find(wireid) != mapOfLinks.end()) {
0226         LogTrace("DTHitAssociator") << "\n" << mapOfLinks[wireid].size() << " Links (DTDigiSimLink):";
0227 
0228         for (vector<DTDigiSimLink>::const_iterator hitIT = mapOfLinks[wireid].begin();
0229              hitIT != mapOfLinks[wireid].end();
0230              ++hitIT) {
0231           LogTrace("DTHitAssociator") << "\t digi number = " << hitIT->number() << ", time = " << hitIT->time()
0232                                       << ", SimTrackId = " << hitIT->SimTrackId();
0233         }
0234       }
0235 
0236       if (mapOfDigi.find(wireid) != mapOfDigi.end()) {
0237         LogTrace("DTHitAssociator") << "\n" << mapOfDigi[wireid].size() << " Digis (DTDigi):";
0238         for (vector<DTDigi>::const_iterator hitIT = mapOfDigi[wireid].begin(); hitIT != mapOfDigi[wireid].end();
0239              ++hitIT) {
0240           LogTrace("DTHitAssociator") << "\t digi number = " << hitIT->number() << ", time = " << hitIT->time();
0241         }
0242       }
0243 
0244       LogTrace("DTHitAssociator") << "\n" << (*mapIT).second.size() << " RecHits (DTRecHit1DPair):";
0245       for (vector<DTRecHit1DPair>::const_iterator vIT = (*mapIT).second.begin(); vIT != (*mapIT).second.end(); ++vIT) {
0246         LogTrace("DTHitAssociator") << "\t digi time = " << vIT->digiTime();
0247       }
0248     }
0249   }
0250 }
0251 // end of constructor
0252 
0253 std::vector<DTHitAssociator::SimHitIdpr> DTHitAssociator::associateHitId(const TrackingRecHit &hit) const {
0254   std::vector<SimHitIdpr> simtrackids;
0255   const TrackingRecHit *hitp = &hit;
0256   const DTRecHit1D *dtrechit = dynamic_cast<const DTRecHit1D *>(hitp);
0257 
0258   if (dtrechit) {
0259     simtrackids = associateDTHitId(dtrechit);
0260   } else {
0261     LogTrace("DTHitAssociator") << "*** WARNING in DTHitAssociator::associateHitId, null dynamic_cast "
0262                                    "!";
0263   }
0264   return simtrackids;
0265 }
0266 
0267 std::vector<DTHitAssociator::SimHitIdpr> DTHitAssociator::associateDTHitId(const DTRecHit1D *dtrechit) const {
0268   std::vector<SimHitIdpr> matched;
0269 
0270   DTWireId wireid = dtrechit->wireId();
0271 
0272   if (config_.associatorByWire) {
0273     // matching based on DTWireId : take only "valid" SimHits on that wire
0274 
0275     auto found = mapOfSimHit.find(wireid);
0276     if (found != mapOfSimHit.end()) {
0277       for (vector<PSimHit_withFlag>::const_iterator hitIT = found->second.begin(); hitIT != found->second.end();
0278            ++hitIT) {
0279         bool valid_hit = hitIT->second;
0280         PSimHit this_hit = hitIT->first;
0281 
0282         if (valid_hit) {
0283           SimHitIdpr currentId(this_hit.trackId(), this_hit.eventId());
0284           matched.push_back(currentId);
0285         }
0286       }
0287     }
0288   }
0289 
0290   else {
0291     // matching based on DTDigiSimLink
0292 
0293     float theTime = dtrechit->digiTime();
0294     int theNumber(-1);
0295 
0296     auto found = mapOfLinks.find(wireid);
0297 
0298     if (found != mapOfLinks.end()) {
0299       // DTDigiSimLink::time() is set equal to DTDigi::time() only for the
0300       // DTDigiSimLink corresponding to the digitizied PSimHit other
0301       // DTDigiSimLinks associated to the same DTDigi are identified by having
0302       // the same DTDigiSimLink::number()
0303 
0304       // first find the associated digi Number
0305       for (vector<DTDigiSimLink>::const_iterator linkIT = found->second.begin(); linkIT != found->second.end();
0306            ++linkIT) {
0307         float digitime = linkIT->time();
0308         if (fabs(digitime - theTime) < 0.1) {
0309           theNumber = linkIT->number();
0310         }
0311       }
0312 
0313       // then get all the DTDigiSimLinks with that digi Number (corresponding to
0314       // valid SimHits
0315       //  within a time window of the order of the time resolution, specified in
0316       //  the DTDigitizer)
0317       for (vector<DTDigiSimLink>::const_iterator linkIT = found->second.begin(); linkIT != found->second.end();
0318            ++linkIT) {
0319         int digiNr = linkIT->number();
0320         if (digiNr == theNumber) {
0321           SimHitIdpr currentId(linkIT->SimTrackId(), linkIT->eventId());
0322           matched.push_back(currentId);
0323         }
0324       }
0325 
0326     } else {
0327       LogTrace("DTHitAssociator") << "*** ERROR in DTHitAssociator::associateDTHitId - DTRecHit1D: " << *dtrechit
0328                                   << " has no associated DTDigiSimLink !" << endl;
0329       return matched;
0330     }
0331   }
0332 
0333   return matched;
0334 }
0335 
0336 std::vector<PSimHit> DTHitAssociator::associateHit(const TrackingRecHit &hit) const {
0337   std::vector<PSimHit> simhits;
0338   std::vector<SimHitIdpr> simtrackids;
0339 
0340   const TrackingRecHit *hitp = &hit;
0341   const DTRecHit1D *dtrechit = dynamic_cast<const DTRecHit1D *>(hitp);
0342 
0343   if (dtrechit) {
0344     DTWireId wireid = dtrechit->wireId();
0345 
0346     if (config_.associatorByWire) {
0347       // matching based on DTWireId : take only "valid" SimHits on that wire
0348 
0349       auto found = mapOfSimHit.find(wireid);
0350       if (found != mapOfSimHit.end()) {
0351         for (vector<PSimHit_withFlag>::const_iterator hitIT = found->second.begin(); hitIT != found->second.end();
0352              ++hitIT) {
0353           bool valid_hit = hitIT->second;
0354           if (valid_hit)
0355             simhits.push_back(hitIT->first);
0356         }
0357       }
0358     } else {
0359       // matching based on DTDigiSimLink
0360 
0361       simtrackids = associateDTHitId(dtrechit);
0362 
0363       for (vector<SimHitIdpr>::const_iterator idIT = simtrackids.begin(); idIT != simtrackids.end(); ++idIT) {
0364         uint32_t trId = idIT->first;
0365         EncodedEventId evId = idIT->second;
0366 
0367         auto found = mapOfSimHit.find(wireid);
0368         if (found != mapOfSimHit.end()) {
0369           for (vector<PSimHit_withFlag>::const_iterator hitIT = found->second.begin(); hitIT != found->second.end();
0370                ++hitIT) {
0371             if (hitIT->first.trackId() == trId && hitIT->first.eventId() == evId)
0372               simhits.push_back(hitIT->first);
0373           }
0374         }
0375       }
0376     }
0377 
0378   } else {
0379     LogTrace("DTHitAssociator") << "*** WARNING in DTHitAssociator::associateHit, null dynamic_cast !";
0380   }
0381   return simhits;
0382 }
0383 
0384 bool DTHitAssociator::SimHitOK(const edm::ESHandle<DTGeometry> &muongeom, const PSimHit &simhit) {
0385   bool result = true;
0386 
0387   DTWireId wireid(simhit.detUnitId());
0388   const DTLayer *dtlayer = muongeom->layer(wireid.layerId());
0389   LocalPoint entryP = simhit.entryPoint();
0390   LocalPoint exitP = simhit.exitPoint();
0391   const DTTopology &topo = dtlayer->specificTopology();
0392   float xwire = topo.wirePosition(wireid.wire());
0393   float xEntry = entryP.x() - xwire;
0394   float xExit = exitP.x() - xwire;
0395   DTTopology::Side entrySide = topo.onWhichBorder(xEntry, entryP.y(), entryP.z());
0396   DTTopology::Side exitSide = topo.onWhichBorder(xExit, exitP.y(), exitP.z());
0397 
0398   bool noParametrisation =
0399       ((entrySide == DTTopology::none || exitSide == DTTopology::none) || (entrySide == exitSide) ||
0400        ((entrySide == DTTopology::xMin && exitSide == DTTopology::xMax) ||
0401         (entrySide == DTTopology::xMax && exitSide == DTTopology::xMin)));
0402 
0403   // discard hits where parametrization can not be used
0404   if (noParametrisation) {
0405     result = false;
0406     return result;
0407   }
0408 
0409   float x;
0410   LocalPoint hitPos = simhit.localPosition();
0411 
0412   if (fabs(hitPos.z()) < 0.002) {
0413     // hit center within 20 um from z=0, no need to extrapolate.
0414     x = hitPos.x() - xwire;
0415   } else {
0416     x = xEntry - (entryP.z() * (xExit - xEntry)) / (exitP.z() - entryP.z());
0417   }
0418 
0419   // discard hits where x is out of range of the parametrization (|x|>2.1 cm)
0420   x *= 10.;  // cm -> mm
0421   if (fabs(x) > 21.)
0422     result = false;
0423 
0424   return result;
0425 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!