Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-06 02:54:13

0001 /*
0002  * StubsSimHitsMatching.cc
0003  *
0004  *  Created on: Jan 13, 2021
0005  *      Author: kbunkow
0006  */
0007 
0008 #include "L1Trigger/L1TMuonOverlapPhase1/interface/Tools/StubsSimHitsMatcher.h"
0009 #include "L1Trigger/L1TMuonOverlapPhase1/interface/Omtf/OmtfName.h"
0010 #include "L1Trigger/L1TMuonOverlapPhase1/interface/Omtf/OMTFinputMaker.h"
0011 #include "L1Trigger/L1TMuonOverlapPhase1/interface/AngleConverterBase.h"
0012 
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 
0016 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0017 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0018 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0019 
0020 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0021 #include "DataFormats/Common/interface/Ptr.h"
0022 
0023 #include "DataFormats/Common/interface/DetSetVector.h"
0024 #include "SimDataFormats/RPCDigiSimLink/interface/RPCDigiSimLink.h"
0025 #include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h"
0026 #include "DataFormats/MuonData/interface/MuonDigiCollection.h"
0027 #include "SimDataFormats/DigiSimLinks/interface/DTDigiSimLink.h"
0028 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0029 
0030 #include "FWCore/ServiceRegistry/interface/Service.h"
0031 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0032 
0033 #include <cmath>
0034 
0035 StubsSimHitsMatcher::StubsSimHitsMatcher(const edm::ParameterSet& edmCfg,
0036                                          const OMTFConfiguration* omtfConfig,
0037                                          const MuonGeometryTokens& muonGeometryTokens)
0038     : omtfConfig(omtfConfig), muonGeometryTokens(muonGeometryTokens) {
0039   rpcSimHitsInputTag = edmCfg.getParameter<edm::InputTag>("rpcSimHitsInputTag");
0040   cscSimHitsInputTag = edmCfg.getParameter<edm::InputTag>("cscSimHitsInputTag");
0041   dtSimHitsInputTag = edmCfg.getParameter<edm::InputTag>("dtSimHitsInputTag");
0042 
0043   rpcDigiSimLinkInputTag = edmCfg.getParameter<edm::InputTag>("rpcDigiSimLinkInputTag");
0044   cscStripDigiSimLinksInputTag = edmCfg.getParameter<edm::InputTag>("cscStripDigiSimLinksInputTag");
0045   dtDigiSimLinksInputTag = edmCfg.getParameter<edm::InputTag>("dtDigiSimLinksInputTag");
0046 
0047   trackingParticleTag = edmCfg.getParameter<edm::InputTag>("trackingParticleTag");
0048 
0049   edm::Service<TFileService> fs;
0050   TFileDirectory subDir = fs->mkdir("StubsSimHitsMatcher");
0051 
0052   allMatchedTracksPdgIds = subDir.make<TH1I>("allMatchedTracksPdgIds", "allMatchedTracksPdgIds", 3, 0, 3);
0053   allMatchedTracksPdgIds->SetCanExtend(TH1::kAllAxes);
0054   bestMatchedTracksPdgIds = subDir.make<TH1I>("bestMatchedTracksPdgIds", "bestMatchedTracksPdgIds", 3, 0, 3);
0055   bestMatchedTracksPdgIds->SetCanExtend(TH1::kAllAxes);
0056 
0057   stubsInLayersCntByPdgId = subDir.make<TH2I>("stubsInLayersCntByPdgId",
0058                                               "stubsInLayersCntByPdgId;;OMTF layer;",
0059                                               3,
0060                                               0,
0061                                               3,
0062                                               omtfConfig->nLayers(),
0063                                               -.5,
0064                                               omtfConfig->nLayers() - 0.5);
0065   stubsInLayersCntByPdgId->SetCanExtend(TH1::kAllAxes);
0066 
0067   firedLayersCntByPdgId = subDir.make<TH2I>("firedLayersCntByPdgId",
0068                                             "firedLayersCntByPdgId",
0069                                             3,
0070                                             0,
0071                                             3,
0072                                             omtfConfig->nLayers(),
0073                                             -.5,
0074                                             omtfConfig->nLayers() - 0.5);
0075   firedLayersCntByPdgId->SetCanExtend(TH1::kAllAxes);
0076 
0077   ptByPdgId = subDir.make<TH2I>("ptByPdgId", "ptByPdgId bestMatched;;pT [GeV];#", 3, 0, 3, 20, 0, 10);
0078   ptByPdgId->SetCanExtend(TH1::kAllAxes);
0079 
0080   rhoByPdgId = subDir.make<TH2I>("rhoByPdgId", "rhoByPdgId bestMatched;;Rho [cm];#", 3, 0, 3, 50, 0, 500);
0081   rhoByPdgId->SetCanExtend(TH1::kAllAxes);
0082 }
0083 
0084 StubsSimHitsMatcher::~StubsSimHitsMatcher() {}
0085 
0086 void StubsSimHitsMatcher::beginRun(edm::EventSetup const& eventSetup) {
0087   if (muonGeometryRecordWatcher.check(eventSetup)) {
0088     _georpc = eventSetup.getHandle(muonGeometryTokens.rpcGeometryEsToken);
0089     _geocsc = eventSetup.getHandle(muonGeometryTokens.cscGeometryEsToken);
0090     _geodt = eventSetup.getHandle(muonGeometryTokens.dtGeometryEsToken);
0091   }
0092 }
0093 
0094 void StubsSimHitsMatcher::match(const edm::Event& iEvent,
0095                                 const l1t::RegionalMuonCand* omtfCand,
0096                                 const AlgoMuonPtr& procMuon,
0097                                 std::ostringstream& ostr) {
0098   ostr << "stubsSimHitsMatching ---------------" << std::endl;
0099 
0100   edm::Handle<edm::PSimHitContainer> rpcSimHitsHandle;
0101   iEvent.getByLabel(rpcSimHitsInputTag, rpcSimHitsHandle);
0102   ostr << "rpcSimHitsHandle: size: " << rpcSimHitsHandle->size() << std::endl;
0103 
0104   edm::Handle<edm::PSimHitContainer> dtSimHitsHandle;
0105   iEvent.getByLabel(dtSimHitsInputTag, dtSimHitsHandle);
0106   ostr << std::endl << "dtSimHitsHandle: size: " << dtSimHitsHandle->size() << std::endl;
0107 
0108   edm::Handle<edm::PSimHitContainer> cscSimHitsHandle;
0109   iEvent.getByLabel(cscSimHitsInputTag, cscSimHitsHandle);
0110   ostr << std::endl << "cscSimHitsHandle: size: " << cscSimHitsHandle->size() << std::endl;
0111 
0112   edm::Handle<edm::DetSetVector<RPCDigiSimLink> > rpcDigiSimLinkHandle;
0113   iEvent.getByLabel(rpcDigiSimLinkInputTag, rpcDigiSimLinkHandle);
0114   ostr << "rpcDigiSimLinkHandle: size: " << rpcDigiSimLinkHandle->size() << std::endl;
0115 
0116   edm::Handle<edm::DetSetVector<StripDigiSimLink> > cscStripDigiSimLinkHandle;
0117   iEvent.getByLabel(cscStripDigiSimLinksInputTag, cscStripDigiSimLinkHandle);
0118   ostr << "cscStripDigiSimLinkHandle: size: " << cscStripDigiSimLinkHandle->size() << std::endl;
0119 
0120   edm::Handle<MuonDigiCollection<DTLayerId, DTDigiSimLink> > dtDigiSimLinkHandle;
0121   iEvent.getByLabel(dtDigiSimLinksInputTag, dtDigiSimLinkHandle);
0122   //ostr<<"dtDigiSimLinkHandle: size: " << dtDigiSimLinkHandle->size()<<std::endl;
0123 
0124   edm::Handle<TrackingParticleCollection> trackingParticleHandle;
0125   iEvent.getByLabel(trackingParticleTag, trackingParticleHandle);
0126 
0127   if (procMuon->isValid() && omtfCand) {
0128     OmtfName board(omtfCand->processor(), omtfCand->trackFinderType());
0129     auto processorPhiZero = OMTFinputMaker::getProcessorPhiZero(omtfConfig, omtfCand->processor());
0130 
0131     std::set<MatchedTrackInfo> matchedTrackInfos;
0132     ostr << board.name() << " " << *procMuon << std::endl;
0133 
0134     auto& gpResult = procMuon->getGpResult();
0135     for (unsigned int iLogicLayer = 0; iLogicLayer < gpResult.getStubResults().size(); ++iLogicLayer) {
0136       auto& stub = gpResult.getStubResults()[iLogicLayer].getMuonStub();
0137       if (stub && gpResult.isLayerFired(iLogicLayer)) {
0138         if (omtfConfig->isBendingLayer(iLogicLayer))
0139           continue;
0140 
0141         DetId stubDetId(stub->detId);
0142         if (stubDetId.det() != DetId::Muon) {
0143           edm::LogError("l1tOmtfEventPrint")
0144               << "!!!!!!!!!!!!!!!!!!!!!!!!  PROBLEM: hit in unknown Det, detID: " << stubDetId.det() << std::endl;
0145           continue;
0146         }
0147 
0148         auto stubGlobalPhi = omtfConfig->procHwPhiToGlobalPhi(stub->phiHw, processorPhiZero);
0149         ostr << (*stub) << "\nstubGlobalPhi " << stubGlobalPhi << std::endl;
0150 
0151         int matchedMuonHits = 0;
0152         int matchedNotMuonHits = 0;
0153 
0154         switch (stubDetId.subdetId()) {
0155           case MuonSubdetId::RPC: {
0156             RPCDetId rpcDetId(stubDetId);
0157 
0158             for (auto& simHit : *(rpcSimHitsHandle.product())) {
0159               if (stubDetId.rawId() == simHit.detUnitId()) {
0160                 const RPCRoll* roll = _georpc->roll(rpcDetId);
0161                 auto strip = roll->strip(simHit.localPosition());
0162                 double simHitStripGlobalPhi = (roll->toGlobal(roll->centreOfStrip((int)strip))).phi();
0163 
0164                 if (abs(stubGlobalPhi - simHitStripGlobalPhi) < 0.02) {
0165                   if (abs(simHit.particleType()) == 13)
0166                     matchedMuonHits++;
0167                   else {
0168                     matchedNotMuonHits++;
0169                   }
0170 
0171                   matchedTrackInfos.insert(MatchedTrackInfo(simHit.eventId().event(), simHit.trackId()));
0172                 }
0173 
0174                 ostr << " simHitStripGlobalPhi " << std::setw(10) << simHitStripGlobalPhi << " strip " << strip
0175                      << " particleType: " << simHit.particleType() << " event: " << simHit.eventId().event()
0176                      << " trackId " << simHit.trackId() << " processType " << simHit.processType() << " detUnitId "
0177                      << simHit.detUnitId() << " "
0178                      << rpcDetId
0179                      //<<" phiAtEntry "<<simHit.phiAtEntry()
0180                      //<<" thetaAtEntry "<<simHit.thetaAtEntry()
0181                      //<<" localPosition: phi "<<simHit.localPosition().phi()<<" eta "<<simHit.localPosition().eta()
0182                      << " entryPoint: x " << std::setw(10) << simHit.entryPoint().x() << " y " << std::setw(10)
0183                      << simHit.entryPoint().y() << " timeOfFlight " << simHit.timeOfFlight() << std::endl;
0184               }
0185             }
0186 
0187             auto rpcDigiSimLinkDetSet = rpcDigiSimLinkHandle.product()->find(stub->detId);
0188 
0189             if (rpcDigiSimLinkDetSet != rpcDigiSimLinkHandle.product()->end()) {
0190               ostr << "rpcDigiSimLinkDetSet: detId " << rpcDigiSimLinkDetSet->detId() << " size "
0191                    << rpcDigiSimLinkDetSet->size() << "\n";
0192               for (auto& rpcDigiSimLink : *rpcDigiSimLinkDetSet) {
0193                 const RPCRoll* roll = _georpc->roll(rpcDetId);
0194                 auto strip = rpcDigiSimLink.getStrip();
0195                 double simHitStripGlobalPhi = (roll->toGlobal(roll->centreOfStrip((int)strip))).phi();
0196 
0197                 if (abs(stubGlobalPhi - simHitStripGlobalPhi) < 0.02) {
0198                   auto matchedTrackInfo = matchedTrackInfos.insert(
0199                       MatchedTrackInfo(rpcDigiSimLink.getEventId().event(), rpcDigiSimLink.getTrackId()));
0200                   matchedTrackInfo.first->matchedDigiCnt.at(iLogicLayer)++;
0201                 }
0202 
0203                 ostr << " simHitStripGlobalPhi " << std::setw(10) << simHitStripGlobalPhi << " strip " << strip
0204                      << " particleType: " << rpcDigiSimLink.getParticleType()
0205                      << " event: " << rpcDigiSimLink.getEventId().event() << " trackId " << rpcDigiSimLink.getTrackId()
0206                      << " processType " << rpcDigiSimLink.getProcessType() << " detUnitId "
0207                      << rpcDigiSimLink.getDetUnitId() << " "
0208                      << rpcDetId
0209                      //<<" phiAtEntry "<<simHit.phiAtEntry()
0210                      //<<" thetaAtEntry "<<simHit.thetaAtEntry()
0211                      //<<" localPosition: phi "<<simHit.localPosition().phi()<<" eta "<<simHit.localPosition().eta()
0212                      << " entryPoint: x " << std::setw(10) << rpcDigiSimLink.getEntryPoint().x() << " y "
0213                      << std::setw(10) << rpcDigiSimLink.getEntryPoint().y() << " timeOfFlight "
0214                      << rpcDigiSimLink.getTimeOfFlight() << std::endl;
0215               }
0216             }
0217 
0218             break;
0219           }  //----------------------------------------------------------------------
0220           case MuonSubdetId::DT: {
0221             //DTChamberId dt(stubDetId);
0222             for (auto& simHit : *(dtSimHitsHandle.product())) {
0223               const DTLayer* layer = _geodt->layer(DTLayerId(simHit.detUnitId()));
0224               const DTChamber* chamber = layer->chamber();
0225               if (stubDetId.rawId() == chamber->id().rawId()) {
0226                 //auto strip = layer->geometry()->strip(simHit.localPosition());
0227                 auto simHitGlobalPoint = layer->toGlobal(simHit.localPosition());
0228 
0229                 ostr << " simHitGlobalPoint.phi " << std::setw(10)
0230                      << simHitGlobalPoint.phi()
0231                      //<<" strip "<<strip
0232                      << " particleType: " << simHit.particleType() << " event: " << simHit.eventId().event()
0233                      << " trackId " << simHit.trackId() << " processType " << simHit.processType() << " detUnitId "
0234                      << simHit.detUnitId() << " "
0235                      << layer->id()
0236                      //<<" phiAtEntry "<<simHit.phiAtEntry()
0237                      //<<" thetaAtEntry "<<simHit.thetaAtEntry()
0238                      //<<" localPosition: phi "<<simHit.localPosition().phi()<<" eta "<<simHit.localPosition().eta()
0239                      << " localPosition: x " << std::setw(10) << simHit.localPosition().x() << " y " << std::setw(10)
0240                      << simHit.localPosition().y() << " timeOfFlight " << simHit.timeOfFlight() << std::endl;
0241               }
0242             }
0243 
0244             auto chamber = _geodt->chamber(DTLayerId(stub->detId));
0245             for (auto superlayer : chamber->superLayers()) {
0246               if (superlayer->id().superLayer() == 2)  //we skip the theta layer
0247                 continue;
0248               for (auto layer : superlayer->layers()) {
0249                 auto dtDigiSimLinks = dtDigiSimLinkHandle.product()->get(layer->id());
0250                 ostr << "dt layer " << layer->id() << "\n";
0251                 auto dtDigiSimLink = dtDigiSimLinks.first;
0252 
0253                 for (; dtDigiSimLink != dtDigiSimLinks.second; dtDigiSimLink++) {
0254                   //const RPCRoll* roll = _georpc->roll(rpcDetId);
0255                   auto wire = dtDigiSimLink->wire();
0256 
0257                   auto wireX = layer->specificTopology().wirePosition(wire);
0258 
0259                   LocalPoint point(wireX, 0, 0);
0260                   auto digiWireGlobal = layer->toGlobal(point);
0261 
0262                   if (abs(stubGlobalPhi - digiWireGlobal.phi()) < 0.03) {
0263                     auto matchedTrackInfo = matchedTrackInfos.insert(
0264                         MatchedTrackInfo(dtDigiSimLink->eventId().event(), dtDigiSimLink->SimTrackId()));
0265                     matchedTrackInfo.first->matchedDigiCnt.at(iLogicLayer)++;
0266                   }
0267 
0268                   ostr
0269                       << " digiWireGlobalPhi " << std::setw(10) << digiWireGlobal.phi() << " wire "
0270                       << wire
0271                       //<<" particleType: "<<cscDigiSimLink.getParticleType()
0272                       << " event: " << dtDigiSimLink->eventId().event() << " trackId "
0273                       << dtDigiSimLink->SimTrackId()
0274                       //<<" CFposition "<<cscDigiSimLink.CFposition() is 0
0275                       //the rest is not available in the StripDigiSimLink, maybe the SimHit must be found in the  CrossingFrame vector(???) with CFposition()
0276                       //<<" processType "<<cscDigiSimLink.getProcessType()
0277                       //<<" detUnitId "<<cscDigiSimLink.getDetUnitId()<<" "<<rpcDetId
0278                       //<<" phiAtEntry "<<simHit.phiAtEntry()
0279                       //<<" thetaAtEntry "<<simHit.thetaAtEntry()
0280                       //<<" localPosition: phi "<<simHit.localPosition().phi()<<" eta "<<simHit.localPosition().eta()
0281                       //<<" entryPoint: x "<<std::setw(10)<<rpcDigiSimLink.getEntryPoint().x()<<" y "<<std::setw(10)<<rpcDigiSimLink.getEntryPoint().y()
0282                       //<<" timeOfFlight "<<rpcDigiSimLink.getTimeOfFlight()
0283                       << std::endl;
0284                 }
0285               }
0286             }
0287 
0288             break;
0289           }  //----------------------------------------------------------------------
0290           case MuonSubdetId::CSC: {
0291             //CSCDetId csc(stubDetId);
0292             for (auto& simHit : *(cscSimHitsHandle.product())) {
0293               const CSCLayer* layer = _geocsc->layer(CSCDetId(simHit.detUnitId()));
0294               auto chamber = layer->chamber();
0295               if (stubDetId.rawId() == chamber->id().rawId()) {
0296                 auto simHitStrip = layer->geometry()->strip(simHit.localPosition());
0297                 auto simHitGlobalPoint = layer->toGlobal(simHit.localPosition());
0298                 auto simHitStripGlobalPhi = layer->centerOfStrip(round(simHitStrip)).phi();
0299 
0300                 ostr << " simHit: gloablPoint phi " << simHitGlobalPoint.phi() << " stripGlobalPhi "
0301                      << simHitStripGlobalPhi.phi() << " strip " << simHitStrip
0302                      << " particleType: " << simHit.particleType() << " event: " << simHit.eventId().event()
0303                      << " trackId " << simHit.trackId() << " processType " << simHit.processType() << " detUnitId "
0304                      << simHit.detUnitId() << " "
0305                      << layer->id()
0306                      //<<" phiAtEntry "<<simHit.phiAtEntry()
0307                      //<<" thetaAtEntry "<<simHit.thetaAtEntry()
0308                      << " timeOfFlight "
0309                      << simHit.timeOfFlight()
0310                      //<<" localPosition: phi "<<simHit.localPosition().phi()<<" eta "<<simHit.localPosition().eta()
0311                      << " x " << simHit.localPosition().x() << " y " << simHit.localPosition().y()
0312 
0313                      << std::endl;
0314               }
0315             }
0316 
0317             auto chamber = _geocsc->chamber(CSCDetId(stub->detId));
0318             for (auto* layer : chamber->layers()) {
0319               auto cscDigiSimLinkDetSet = cscStripDigiSimLinkHandle.product()->find(layer->id());
0320               if (cscDigiSimLinkDetSet != cscStripDigiSimLinkHandle.product()->end()) {
0321                 ostr << "cscDigiSimLinkDetSet: detId " << cscDigiSimLinkDetSet->detId() << " " << layer->id()
0322                      << " size " << cscDigiSimLinkDetSet->size() << "\n";
0323                 for (auto& cscDigiSimLink : *cscDigiSimLinkDetSet) {
0324                   //const RPCRoll* roll = _georpc->roll(rpcDetId);
0325                   auto strip = cscDigiSimLink.channel();
0326                   auto digiStripGlobalPhi = layer->centerOfStrip(strip).phi();
0327 
0328                   if (abs(stubGlobalPhi - digiStripGlobalPhi) < 0.03) {
0329                     auto matchedTrackInfo = matchedTrackInfos.insert(
0330                         MatchedTrackInfo(cscDigiSimLink.eventId().event(), cscDigiSimLink.SimTrackId()));
0331                     matchedTrackInfo.first->matchedDigiCnt.at(iLogicLayer)++;
0332                   }
0333                   ostr
0334                       << " digiStripGlobalPhi " << std::setw(10) << digiStripGlobalPhi << " strip "
0335                       << strip
0336                       //<<" particleType: "<<cscDigiSimLink.getParticleType()
0337                       << " event: " << cscDigiSimLink.eventId().event() << " trackId "
0338                       << cscDigiSimLink.SimTrackId()
0339                       //<<" CFposition "<<cscDigiSimLink.CFposition() is 0
0340                       //the rest is not available in the StripDigiSimLink, maybe the SimHit must be found in the  CrossingFrame vector(???) with CFposition()
0341                       //<<" processType "<<cscDigiSimLink.getProcessType()
0342                       //<<" detUnitId "<<cscDigiSimLink.getDetUnitId()<<" "<<rpcDetId
0343                       //<<" phiAtEntry "<<simHit.phiAtEntry()
0344                       //<<" thetaAtEntry "<<simHit.thetaAtEntry()
0345                       //<<" localPosition: phi "<<simHit.localPosition().phi()<<" eta "<<simHit.localPosition().eta()
0346                       //<<" entryPoint: x "<<std::setw(10)<<rpcDigiSimLink.getEntryPoint().x()<<" y "<<std::setw(10)<<rpcDigiSimLink.getEntryPoint().y()
0347                       //<<" timeOfFlight "<<rpcDigiSimLink.getTimeOfFlight()
0348                       << std::endl;
0349                 }
0350               } else {
0351                 ostr << "cscDigiSimLinkDetSet not found for detId " << layer->id();
0352               }
0353             }
0354 
0355             break;
0356           }  //end of CSC case
0357         }    //end of switch
0358         ostr << "" << std::endl;
0359       }
0360     }
0361 
0362     ostr << board.name() << " " << *procMuon << std::endl;
0363     ostr << procMuon->getGpResult() << std::endl << std::endl;
0364 
0365     int maxMatchedStubs = 0;
0366     const TrackingParticle* bestMatchedPart = nullptr;
0367     for (auto matchedTrackInfo : matchedTrackInfos) {
0368       ostr << "matchedTrackInfo eventNum " << matchedTrackInfo.eventNum << " trackId " << matchedTrackInfo.trackId
0369            << "\n";
0370 
0371       const TrackingParticle* matchedPart = nullptr;
0372       for (auto& trackingParticle :
0373            *trackingParticleHandle.product()) {  //finding the trackingParticle corresponding to the matchedTrackInfo
0374         if (trackingParticle.eventId().event() == matchedTrackInfo.eventNum &&
0375             trackingParticle.g4Tracks().at(0).trackId() == matchedTrackInfo.trackId) {
0376           allMatchedTracksPdgIds->Fill(to_string(trackingParticle.pdgId()).c_str(), 1);
0377           matchedPart = &trackingParticle;
0378 
0379           ostr << "trackingParticle: pdgId " << std::setw(3) << trackingParticle.pdgId() << " event " << std::setw(4)
0380                << trackingParticle.eventId().event() << " trackId " << std::setw(8)
0381                << trackingParticle.g4Tracks().at(0).trackId() << " pt " << std::setw(9) << trackingParticle.pt()
0382                << " eta " << std::setw(9) << trackingParticle.momentum().eta() << " phi " << std::setw(9)
0383                << trackingParticle.momentum().phi() << std::endl;
0384 
0385           if (trackingParticle.parentVertex().isNonnull()) {
0386             ostr << "parentVertex Rho " << trackingParticle.parentVertex()->position().Rho() << " z "
0387                  << trackingParticle.parentVertex()->position().z() << " R "
0388                  << trackingParticle.parentVertex()->position().R() << " eta "
0389                  << trackingParticle.parentVertex()->position().eta() << " phi "
0390                  << trackingParticle.parentVertex()->position().phi() << std::endl;
0391 
0392             for (auto& parentTrack : trackingParticle.parentVertex()->sourceTracks()) {
0393               ostr << "parentTrack:      pdgId " << std::setw(3) << parentTrack->pdgId() << " event " << std::setw(4)
0394                    << parentTrack->eventId().event() << " trackId " << std::setw(8)
0395                    << parentTrack->g4Tracks().at(0).trackId() << " pt " << std::setw(9) << parentTrack->pt() << " eta "
0396                    << std::setw(9) << parentTrack->momentum().eta() << " phi " << std::setw(9)
0397                    << parentTrack->momentum().phi() << std::endl;
0398             }
0399           }
0400 
0401           break;
0402         }
0403       }
0404 
0405       int matchedStubsCnt = 0;
0406       for (unsigned int iLayer = 0; iLayer < matchedTrackInfo.matchedDigiCnt.size(); iLayer++) {
0407         ostr << " matchedDigiCnt in layer " << iLayer << " " << matchedTrackInfo.matchedDigiCnt[iLayer] << "\n";
0408         if (matchedTrackInfo.matchedDigiCnt[iLayer] > 0) {
0409           matchedStubsCnt++;
0410           if (matchedPart) {
0411             string str = to_string(matchedPart->pdgId());
0412             stubsInLayersCntByPdgId->Fill(str.c_str(), iLayer, 1);
0413           }
0414         }
0415       }
0416       ostr << std::endl;
0417 
0418       if (maxMatchedStubs < matchedStubsCnt) {
0419         maxMatchedStubs = matchedStubsCnt;
0420         bestMatchedPart = matchedPart;
0421       }
0422     }
0423     if (bestMatchedPart) {
0424       string str = to_string(bestMatchedPart->pdgId());
0425 
0426       bestMatchedTracksPdgIds->Fill(str.c_str(), 1);
0427       firedLayersCntByPdgId->Fill(str.c_str(), gpResult.getFiredLayerCnt(), 1);
0428       ptByPdgId->Fill(str.c_str(), bestMatchedPart->pt(), 1);
0429 
0430       if (bestMatchedPart->parentVertex().isNonnull()) {
0431         rhoByPdgId->Fill(str.c_str(), bestMatchedPart->parentVertex()->position().Rho(), 1);
0432       }
0433 
0434       ostr << "bestMatchedPart: pdgId " << std::setw(3) << bestMatchedPart->pdgId() << " event " << std::setw(4)
0435            << bestMatchedPart->eventId().event() << " trackId " << std::setw(8)
0436            << bestMatchedPart->g4Tracks().at(0).trackId() << "\n\n";
0437     }
0438   }
0439 }
0440 
0441 void StubsSimHitsMatcher::endJob() {
0442   /*edm::LogVerbatim("l1tOmtfEventPrint") <<"allMatchedTracksPdgIds "<<std::endl;
0443   for(auto pdgId : allMatchedTracksPdgIds) {
0444     edm::LogVerbatim("l1tOmtfEventPrint") <<"pdgId "<<std::setw(8)<<pdgId.first<<" muon candidates "<<pdgId.second<<std::endl;
0445   }
0446 
0447   edm::LogVerbatim("l1tOmtfEventPrint") <<"bestMatchedTracksPdgIds "<<std::endl;
0448   for(auto pdgId : bestMatchedTracksPdgIds) {
0449     edm::LogVerbatim("l1tOmtfEventPrint") <<"pdgId "<<std::setw(8) <<pdgId.first<<" muon candidates "<<pdgId.second<<std::endl;
0450   }*/
0451 }