src/Tools/StubsSimHitsMatching.cc

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(), omtfConfig);
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->getGpResultConstr();
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         switch (stubDetId.subdetId()) {
0152           case MuonSubdetId::RPC: {
0153             RPCDetId rpcDetId(stubDetId);
0154
0155             for (auto& simHit : *(rpcSimHitsHandle.product())) {
0156               if (stubDetId.rawId() == simHit.detUnitId()) {
0157                 const RPCRoll* roll = _georpc->roll(rpcDetId);
0158                 auto strip = roll->strip(simHit.localPosition());
0159                 double simHitStripGlobalPhi = (roll->toGlobal(roll->centreOfStrip((int)strip))).phi();
0160
0161                 if (std::abs(stubGlobalPhi - simHitStripGlobalPhi) < 0.02) {
0162                   matchedTrackInfos.insert(MatchedTrackInfo(simHit.eventId().event(), simHit.trackId()));
0163                 }
0164
0165                 ostr << " simHitStripGlobalPhi " << std::setw(10) << simHitStripGlobalPhi << " strip " << strip
0166                      << " particleType: " << simHit.particleType() << " event: " << simHit.eventId().event()
0167                      << " trackId " << simHit.trackId() << " processType " << simHit.processType() << " detUnitId "
0168                      << simHit.detUnitId() << " "
0169                      << rpcDetId
0170                      //<<" phiAtEntry "<<simHit.phiAtEntry()
0171                      //<<" thetaAtEntry "<<simHit.thetaAtEntry()
0172                      //<<" localPosition: phi "<<simHit.localPosition().phi()<<" eta "<<simHit.localPosition().eta()
0173                      << " entryPoint: x " << std::setw(10) << simHit.entryPoint().x() << " y " << std::setw(10)
0174                      << simHit.entryPoint().y() << " timeOfFlight " << simHit.timeOfFlight() << std::endl;
0175               }
0176             }
0177
0178             auto rpcDigiSimLinkDetSet = rpcDigiSimLinkHandle.product()->find(stub->detId);
0179
0180             if (rpcDigiSimLinkDetSet != rpcDigiSimLinkHandle.product()->end()) {
0181               ostr << "rpcDigiSimLinkDetSet: detId " << rpcDigiSimLinkDetSet->detId() << " size "
0182                    << rpcDigiSimLinkDetSet->size() << "\n";
0183               for (auto& rpcDigiSimLink : *rpcDigiSimLinkDetSet) {
0184                 const RPCRoll* roll = _georpc->roll(rpcDetId);
0185                 auto strip = rpcDigiSimLink.getStrip();
0186                 double simHitStripGlobalPhi = (roll->toGlobal(roll->centreOfStrip((int)strip))).phi();
0187
0188                 if (std::abs(stubGlobalPhi - simHitStripGlobalPhi) < 0.02) {
0189                   auto matchedTrackInfo = matchedTrackInfos.insert(
0190                       MatchedTrackInfo(rpcDigiSimLink.getEventId().event(), rpcDigiSimLink.getTrackId()));
0191                   matchedTrackInfo.first->matchedDigiCnt.at(iLogicLayer)++;
0192                 }
0193
0194                 ostr << " simHitStripGlobalPhi " << std::setw(10) << simHitStripGlobalPhi << " strip " << strip
0195                      << " particleType: " << rpcDigiSimLink.getParticleType()
0196                      << " event: " << rpcDigiSimLink.getEventId().event() << " trackId " << rpcDigiSimLink.getTrackId()
0197                      << " processType " << rpcDigiSimLink.getProcessType() << " detUnitId "
0198                      << rpcDigiSimLink.getDetUnitId() << " "
0199                      << rpcDetId
0200                      //<<" phiAtEntry "<<simHit.phiAtEntry()
0201                      //<<" thetaAtEntry "<<simHit.thetaAtEntry()
0202                      //<<" localPosition: phi "<<simHit.localPosition().phi()<<" eta "<<simHit.localPosition().eta()
0203                      << " entryPoint: x " << std::setw(10) << rpcDigiSimLink.getEntryPoint().x() << " y "
0204                      << std::setw(10) << rpcDigiSimLink.getEntryPoint().y() << " timeOfFlight "
0205                      << rpcDigiSimLink.getTimeOfFlight() << std::endl;
0206               }
0207             }
0208
0209             break;
0210           }  //----------------------------------------------------------------------
0211           case MuonSubdetId::DT: {
0212             //DTChamberId dt(stubDetId);
0213             for (auto& simHit : *(dtSimHitsHandle.product())) {
0214               const DTLayer* layer = _geodt->layer(DTLayerId(simHit.detUnitId()));
0215               const DTChamber* chamber = layer->chamber();
0216               if (stubDetId.rawId() == chamber->id().rawId()) {
0217                 //auto strip = layer->geometry()->strip(simHit.localPosition());
0218                 auto simHitGlobalPoint = layer->toGlobal(simHit.localPosition());
0219
0220                 ostr << " simHitGlobalPoint.phi " << std::setw(10)
0221                      << simHitGlobalPoint.phi()
0222                      //<<" strip "<<strip
0223                      << " particleType: " << simHit.particleType() << " event: " << simHit.eventId().event()
0224                      << " trackId " << simHit.trackId() << " processType " << simHit.processType() << " detUnitId "
0225                      << simHit.detUnitId() << " "
0226                      << layer->id()
0227                      //<<" phiAtEntry "<<simHit.phiAtEntry()
0228                      //<<" thetaAtEntry "<<simHit.thetaAtEntry()
0229                      //<<" localPosition: phi "<<simHit.localPosition().phi()<<" eta "<<simHit.localPosition().eta()
0230                      << " localPosition: x " << std::setw(10) << simHit.localPosition().x() << " y " << std::setw(10)
0231                      << simHit.localPosition().y() << " timeOfFlight " << simHit.timeOfFlight() << std::endl;
0232               }
0233             }
0234
0235             auto chamber = _geodt->chamber(DTLayerId(stub->detId));
0236             for (auto superlayer : chamber->superLayers()) {
0237               if (superlayer->id().superLayer() == 2)  //we skip the theta layer
0238                 continue;
0239               for (auto layer : superlayer->layers()) {
0240                 auto dtDigiSimLinks = dtDigiSimLinkHandle.product()->get(layer->id());
0241                 ostr << "dt layer " << layer->id() << "\n";
0242                 auto dtDigiSimLink = dtDigiSimLinks.first;
0243
0244                 for (; dtDigiSimLink != dtDigiSimLinks.second; dtDigiSimLink++) {
0245                   //const RPCRoll* roll = _georpc->roll(rpcDetId);
0246                   auto wire = dtDigiSimLink->wire();
0247
0248                   auto wireX = layer->specificTopology().wirePosition(wire);
0249
0250                   LocalPoint point(wireX, 0, 0);
0251                   auto digiWireGlobal = layer->toGlobal(point);
0252
0253                   if (std::abs(stubGlobalPhi - digiWireGlobal.phi()) < 0.03) {
0254                     auto matchedTrackInfo = matchedTrackInfos.insert(
0255                         MatchedTrackInfo(dtDigiSimLink->eventId().event(), dtDigiSimLink->SimTrackId()));
0256                     matchedTrackInfo.first->matchedDigiCnt.at(iLogicLayer)++;
0257                   }
0258
0259                   ostr
0260                       << " digiWireGlobalPhi " << std::setw(10) << digiWireGlobal.phi() << " wire "
0261                       << wire
0262                       //<<" particleType: "<<cscDigiSimLink.getParticleType()
0263                       << " event: " << dtDigiSimLink->eventId().event() << " trackId "
0264                       << dtDigiSimLink->SimTrackId()
0265                       //<<" CFposition "<<cscDigiSimLink.CFposition() is 0
0266                       //the rest is not available in the StripDigiSimLink, maybe the SimHit must be found in the  CrossingFrame vector(???) with CFposition()
0267                       //<<" processType "<<cscDigiSimLink.getProcessType()
0268                       //<<" detUnitId "<<cscDigiSimLink.getDetUnitId()<<" "<<rpcDetId
0269                       //<<" phiAtEntry "<<simHit.phiAtEntry()
0270                       //<<" thetaAtEntry "<<simHit.thetaAtEntry()
0271                       //<<" localPosition: phi "<<simHit.localPosition().phi()<<" eta "<<simHit.localPosition().eta()
0272                       //<<" entryPoint: x "<<std::setw(10)<<rpcDigiSimLink.getEntryPoint().x()<<" y "<<std::setw(10)<<rpcDigiSimLink.getEntryPoint().y()
0273                       //<<" timeOfFlight "<<rpcDigiSimLink.getTimeOfFlight()
0274                       << std::endl;
0275                 }
0276               }
0277             }
0278
0279             break;
0280           }  //----------------------------------------------------------------------
0281           case MuonSubdetId::CSC: {
0282             //CSCDetId csc(stubDetId);
0283             for (auto& simHit : *(cscSimHitsHandle.product())) {
0284               const CSCLayer* layer = _geocsc->layer(CSCDetId(simHit.detUnitId()));
0285               auto chamber = layer->chamber();
0286               if (stubDetId.rawId() == chamber->id().rawId()) {
0287                 auto simHitStrip = layer->geometry()->strip(simHit.localPosition());
0288                 auto simHitGlobalPoint = layer->toGlobal(simHit.localPosition());
0289                 auto simHitStripGlobalPhi = layer->centerOfStrip(round(simHitStrip)).phi();
0290
0291                 ostr << " simHit: gloablPoint phi " << simHitGlobalPoint.phi() << " stripGlobalPhi "
0292                      << simHitStripGlobalPhi.phi() << " strip " << simHitStrip
0293                      << " particleType: " << simHit.particleType() << " event: " << simHit.eventId().event()
0294                      << " trackId " << simHit.trackId() << " processType " << simHit.processType() << " detUnitId "
0295                      << simHit.detUnitId() << " "
0296                      << layer->id()
0297                      //<<" phiAtEntry "<<simHit.phiAtEntry()
0298                      //<<" thetaAtEntry "<<simHit.thetaAtEntry()
0299                      << " timeOfFlight "
0300                      << simHit.timeOfFlight()
0301                      //<<" localPosition: phi "<<simHit.localPosition().phi()<<" eta "<<simHit.localPosition().eta()
0302                      << " x " << simHit.localPosition().x() << " y " << simHit.localPosition().y()
0303
0304                      << std::endl;
0305               }
0306             }
0307
0308             auto chamber = _geocsc->chamber(CSCDetId(stub->detId));
0309             for (auto* layer : chamber->layers()) {
0310               auto cscDigiSimLinkDetSet = cscStripDigiSimLinkHandle.product()->find(layer->id());
0311               if (cscDigiSimLinkDetSet != cscStripDigiSimLinkHandle.product()->end()) {
0312                 ostr << "cscDigiSimLinkDetSet: detId " << cscDigiSimLinkDetSet->detId() << " " << layer->id()
0313                      << " size " << cscDigiSimLinkDetSet->size() << "\n";
0314                 for (auto& cscDigiSimLink : *cscDigiSimLinkDetSet) {
0315                   //const RPCRoll* roll = _georpc->roll(rpcDetId);
0316                   auto strip = cscDigiSimLink.channel();
0317                   auto digiStripGlobalPhi = layer->centerOfStrip(strip).phi();
0318
0319                   if (std::abs(stubGlobalPhi - digiStripGlobalPhi) < 0.03) {
0320                     auto matchedTrackInfo = matchedTrackInfos.insert(
0321                         MatchedTrackInfo(cscDigiSimLink.eventId().event(), cscDigiSimLink.SimTrackId()));
0322                     matchedTrackInfo.first->matchedDigiCnt.at(iLogicLayer)++;
0323                   }
0324                   ostr
0325                       << " digiStripGlobalPhi " << std::setw(10) << digiStripGlobalPhi << " strip "
0326                       << strip
0327                       //<<" particleType: "<<cscDigiSimLink.getParticleType()
0328                       << " event: " << cscDigiSimLink.eventId().event() << " trackId "
0329                       << cscDigiSimLink.SimTrackId()
0330                       //<<" CFposition "<<cscDigiSimLink.CFposition() is 0
0331                       //the rest is not available in the StripDigiSimLink, maybe the SimHit must be found in the  CrossingFrame vector(???) with CFposition()
0332                       //<<" processType "<<cscDigiSimLink.getProcessType()
0333                       //<<" detUnitId "<<cscDigiSimLink.getDetUnitId()<<" "<<rpcDetId
0334                       //<<" phiAtEntry "<<simHit.phiAtEntry()
0335                       //<<" thetaAtEntry "<<simHit.thetaAtEntry()
0336                       //<<" localPosition: phi "<<simHit.localPosition().phi()<<" eta "<<simHit.localPosition().eta()
0337                       //<<" entryPoint: x "<<std::setw(10)<<rpcDigiSimLink.getEntryPoint().x()<<" y "<<std::setw(10)<<rpcDigiSimLink.getEntryPoint().y()
0338                       //<<" timeOfFlight "<<rpcDigiSimLink.getTimeOfFlight()
0339                       << std::endl;
0340                 }
0341               } else {
0342                 ostr << "cscDigiSimLinkDetSet not found for detId " << layer->id();
0343               }
0344             }
0345
0346             break;
0347           }  //end of CSC case
0348         }    //end of switch
0349         ostr << "" << std::endl;
0350       }
0351     }
0352
0353     ostr << board.name() << " " << *procMuon << std::endl;
0354     ostr << procMuon->getGpResultConstr() << std::endl << std::endl;
0355
0356     int maxMatchedStubs = 0;
0357     const TrackingParticle* bestMatchedPart = nullptr;
0358     for (auto matchedTrackInfo : matchedTrackInfos) {
0359       ostr << "matchedTrackInfo eventNum " << matchedTrackInfo.eventNum << " trackId " << matchedTrackInfo.trackId
0360            << "\n";
0361
0362       const TrackingParticle* matchedPart = nullptr;
0363       for (auto& trackingParticle :
0364            *trackingParticleHandle.product()) {  //finding the trackingParticle corresponding to the matchedTrackInfo
0365         if (trackingParticle.eventId().event() == matchedTrackInfo.eventNum &&
0366             trackingParticle.g4Tracks().at(0).trackId() == matchedTrackInfo.trackId) {
0367           allMatchedTracksPdgIds->Fill(to_string(trackingParticle.pdgId()).c_str(), 1);
0368           matchedPart = &trackingParticle;
0369
0370           ostr << "trackingParticle: pdgId " << std::setw(3) << trackingParticle.pdgId() << " event " << std::setw(4)
0371                << trackingParticle.eventId().event() << " trackId " << std::setw(8)
0372                << trackingParticle.g4Tracks().at(0).trackId() << " pt " << std::setw(9) << trackingParticle.pt()
0373                << " eta " << std::setw(9) << trackingParticle.momentum().eta() << " phi " << std::setw(9)
0374                << trackingParticle.momentum().phi() << std::endl;
0375
0376           if (trackingParticle.parentVertex().isNonnull()) {
0377             ostr << "parentVertex Rho " << trackingParticle.parentVertex()->position().Rho() << " z "
0378                  << trackingParticle.parentVertex()->position().z() << " R "
0379                  << trackingParticle.parentVertex()->position().R() << " eta "
0380                  << trackingParticle.parentVertex()->position().eta() << " phi "
0381                  << trackingParticle.parentVertex()->position().phi() << std::endl;
0382
0383             for (auto& parentTrack : trackingParticle.parentVertex()->sourceTracks()) {
0384               ostr << "parentTrack:      pdgId " << std::setw(3) << parentTrack->pdgId() << " event " << std::setw(4)
0385                    << parentTrack->eventId().event() << " trackId " << std::setw(8)
0386                    << parentTrack->g4Tracks().at(0).trackId() << " pt " << std::setw(9) << parentTrack->pt() << " eta "
0387                    << std::setw(9) << parentTrack->momentum().eta() << " phi " << std::setw(9)
0388                    << parentTrack->momentum().phi() << std::endl;
0389             }
0390           }
0391
0392           break;
0393         }
0394       }
0395
0396       int matchedStubsCnt = 0;
0397       for (unsigned int iLayer = 0; iLayer < matchedTrackInfo.matchedDigiCnt.size(); iLayer++) {
0398         ostr << " matchedDigiCnt in layer " << iLayer << " " << matchedTrackInfo.matchedDigiCnt[iLayer] << "\n";
0399         if (matchedTrackInfo.matchedDigiCnt[iLayer] > 0) {
0400           matchedStubsCnt++;
0401           if (matchedPart) {
0402             string str = to_string(matchedPart->pdgId());
0403             stubsInLayersCntByPdgId->Fill(str.c_str(), iLayer, 1);
0404           }
0405         }
0406       }
0407       ostr << std::endl;
0408
0409       if (maxMatchedStubs < matchedStubsCnt) {
0410         maxMatchedStubs = matchedStubsCnt;
0411         bestMatchedPart = matchedPart;
0412       }
0413     }
0414     if (bestMatchedPart) {
0415       string str = to_string(bestMatchedPart->pdgId());
0416
0417       bestMatchedTracksPdgIds->Fill(str.c_str(), 1);
0418       firedLayersCntByPdgId->Fill(str.c_str(), gpResult.getFiredLayerCnt(), 1);
0419       ptByPdgId->Fill(str.c_str(), bestMatchedPart->pt(), 1);
0420
0421       if (bestMatchedPart->parentVertex().isNonnull()) {
0422         rhoByPdgId->Fill(str.c_str(), bestMatchedPart->parentVertex()->position().Rho(), 1);
0423       }
0424
0425       ostr << "bestMatchedPart: pdgId " << std::setw(3) << bestMatchedPart->pdgId() << " event " << std::setw(4)
0426            << bestMatchedPart->eventId().event() << " trackId " << std::setw(8)
0427            << bestMatchedPart->g4Tracks().at(0).trackId() << "\n\n";
0428     }
0429   }
0430 }
0431
0432 void StubsSimHitsMatcher::endJob() {
0433   /*edm::LogVerbatim("l1tOmtfEventPrint") <<"allMatchedTracksPdgIds "<<std::endl;
0434   for(auto pdgId : allMatchedTracksPdgIds) {
0435     edm::LogVerbatim("l1tOmtfEventPrint") <<"pdgId "<<std::setw(8)<<pdgId.first<<" muon candidates "<<pdgId.second<<std::endl;
0436   }
0437
0438   edm::LogVerbatim("l1tOmtfEventPrint") <<"bestMatchedTracksPdgIds "<<std::endl;
0439   for(auto pdgId : bestMatchedTracksPdgIds) {
0440     edm::LogVerbatim("l1tOmtfEventPrint") <<"pdgId "<<std::setw(8) <<pdgId.first<<" muon candidates "<<pdgId.second<<std::endl;
0441   }*/
0442 }