Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:33:21 

0001 #include "Validation/RecoTrack/interface/MultiTrackValidator.h"
0002 #include "Validation/RecoTrack/interface/trackFromSeedFitFailed.h"
0003 
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/Utilities/interface/transform.h"
0007 
0008 #include "DataFormats/TrackReco/interface/Track.h"
0009 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0010 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0011 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0012 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0013 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0014 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0015 #include "SimTracker/TrackAssociation/interface/CosmicParametersDefinerForTP.h"
0016 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0017 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0018 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0019 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0020 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0021 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
0022 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0023 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0024 #include "SimTracker/TrackAssociation/interface/TrackingParticleIP.h"
0025 
0026 #include "DataFormats/TrackReco/interface/DeDxData.h"
0027 #include "DataFormats/Common/interface/ValueMap.h"
0028 #include "DataFormats/Common/interface/Ref.h"
0029 #include "CommonTools/Utils/interface/associationMapFilterValues.h"
0030 #include "FWCore/Utilities/interface/IndexSet.h"
0031 #include <type_traits>
0032 
0033 #include "TMath.h"
0034 #include <TF1.h>
0035 #include "DataFormats/Math/interface/deltaR.h"
0036 #include "DataFormats/Math/interface/PtEtaPhiMass.h"
0037 //#include <iostream>
0038 
0039 using namespace std;
0040 using namespace edm;
0041 
0042 typedef edm::Ref<edm::HepMCProduct, HepMC::GenParticle> GenParticleRef;
0043 namespace {
0044   bool trackSelected(unsigned char mask, unsigned char qual) { return mask & 1 << qual; }
0045 
0046 }  // namespace
0047 
0048 MultiTrackValidator::MultiTrackValidator(const edm::ParameterSet& pset)
0049     : tTopoEsToken(esConsumes()),
0050       parametersDefinerIsCosmic_(pset.getParameter<std::string>("parametersDefiner") == "CosmicParametersDefinerForTP"),
0051       associators(pset.getUntrackedParameter<std::vector<edm::InputTag>>("associators")),
0052       label(pset.getParameter<std::vector<edm::InputTag>>("label")),
0053       ignoremissingtkcollection_(pset.getUntrackedParameter<bool>("ignoremissingtrackcollection", false)),
0054       useAssociators_(pset.getParameter<bool>("UseAssociators")),
0055       calculateDrSingleCollection_(pset.getUntrackedParameter<bool>("calculateDrSingleCollection")),
0056       doPlotsOnlyForTruePV_(pset.getUntrackedParameter<bool>("doPlotsOnlyForTruePV")),
0057       doSummaryPlots_(pset.getUntrackedParameter<bool>("doSummaryPlots")),
0058       doSimPlots_(pset.getUntrackedParameter<bool>("doSimPlots")),
0059       doSimTrackPlots_(pset.getUntrackedParameter<bool>("doSimTrackPlots")),
0060       doRecoTrackPlots_(pset.getUntrackedParameter<bool>("doRecoTrackPlots")),
0061       dodEdxPlots_(pset.getUntrackedParameter<bool>("dodEdxPlots")),
0062       doPVAssociationPlots_(pset.getUntrackedParameter<bool>("doPVAssociationPlots")),
0063       doSeedPlots_(pset.getUntrackedParameter<bool>("doSeedPlots")),
0064       doMVAPlots_(pset.getUntrackedParameter<bool>("doMVAPlots")),
0065       applyTPSelToSimMatch_(pset.getParameter<bool>("applyTPSelToSimMatch")),
0066       simPVMaxZ_(pset.getUntrackedParameter<double>("simPVMaxZ")) {
0067   if (not(pset.getParameter<edm::InputTag>("cores").label().empty())) {
0068     cores_ = consumes<edm::View<reco::Candidate>>(pset.getParameter<edm::InputTag>("cores"));
0069   }
0070   if (label.empty()) {
0071     // Disable prefetching of everything if there are no track collections
0072     return;
0073   }
0074 
0075   const edm::InputTag& label_tp_effic_tag = pset.getParameter<edm::InputTag>("label_tp_effic");
0076   const edm::InputTag& label_tp_fake_tag = pset.getParameter<edm::InputTag>("label_tp_fake");
0077 
0078   if (pset.getParameter<bool>("label_tp_effic_refvector")) {
0079     label_tp_effic_refvector = consumes<TrackingParticleRefVector>(label_tp_effic_tag);
0080   } else {
0081     label_tp_effic = consumes<TrackingParticleCollection>(label_tp_effic_tag);
0082   }
0083   if (pset.getParameter<bool>("label_tp_fake_refvector")) {
0084     label_tp_fake_refvector = consumes<TrackingParticleRefVector>(label_tp_fake_tag);
0085   } else {
0086     label_tp_fake = consumes<TrackingParticleCollection>(label_tp_fake_tag);
0087   }
0088   label_pileupinfo = consumes<std::vector<PileupSummaryInfo>>(pset.getParameter<edm::InputTag>("label_pileupinfo"));
0089   for (const auto& tag : pset.getParameter<std::vector<edm::InputTag>>("sim")) {
0090     simHitTokens_.push_back(consumes<std::vector<PSimHit>>(tag));
0091   }
0092 
0093   std::vector<edm::InputTag> doResolutionPlotsForLabels =
0094       pset.getParameter<std::vector<edm::InputTag>>("doResolutionPlotsForLabels");
0095   doResolutionPlots_.reserve(label.size());
0096   for (auto& itag : label) {
0097     labelToken.push_back(consumes<edm::View<reco::Track>>(itag));
0098     const bool doResol = doResolutionPlotsForLabels.empty() ||
0099                          (std::find(cbegin(doResolutionPlotsForLabels), cend(doResolutionPlotsForLabels), itag) !=
0100                           cend(doResolutionPlotsForLabels));
0101     doResolutionPlots_.push_back(doResol);
0102   }
0103   {  // check for duplicates
0104     auto labelTmp = edm::vector_transform(label, [&](const edm::InputTag& tag) { return tag.label(); });
0105     std::sort(begin(labelTmp), end(labelTmp));
0106     std::string empty;
0107     const std::string* prev = &empty;
0108     for (const std::string& l : labelTmp) {
0109       if (l == *prev) {
0110         throw cms::Exception("Configuration") << "Duplicate InputTag in labels: " << l;
0111       }
0112       prev = &l;
0113     }
0114   }
0115 
0116   edm::InputTag beamSpotTag = pset.getParameter<edm::InputTag>("beamSpot");
0117   bsSrc = consumes<reco::BeamSpot>(beamSpotTag);
0118 
0119   if (parametersDefinerIsCosmic_) {
0120     parametersDefinerTP_ = std::make_unique<CosmicParametersDefinerForTP>(consumesCollector());
0121   } else {
0122     parametersDefinerTP_ = std::make_unique<ParametersDefinerForTP>(beamSpotTag, consumesCollector());
0123   }
0124 
0125   ParameterSet psetForHistoProducerAlgo = pset.getParameter<ParameterSet>("histoProducerAlgoBlock");
0126   histoProducerAlgo_ = std::make_unique<MTVHistoProducerAlgoForTracker>(psetForHistoProducerAlgo, doSeedPlots_);
0127 
0128   dirName_ = pset.getParameter<std::string>("dirName");
0129 
0130   tpNLayersToken_ = consumes<edm::ValueMap<unsigned int>>(pset.getParameter<edm::InputTag>("label_tp_nlayers"));
0131   tpNPixelLayersToken_ =
0132       consumes<edm::ValueMap<unsigned int>>(pset.getParameter<edm::InputTag>("label_tp_npixellayers"));
0133   tpNStripStereoLayersToken_ =
0134       consumes<edm::ValueMap<unsigned int>>(pset.getParameter<edm::InputTag>("label_tp_nstripstereolayers"));
0135 
0136   if (dodEdxPlots_) {
0137     m_dEdx1Tag = consumes<edm::ValueMap<reco::DeDxData>>(pset.getParameter<edm::InputTag>("dEdx1Tag"));
0138     m_dEdx2Tag = consumes<edm::ValueMap<reco::DeDxData>>(pset.getParameter<edm::InputTag>("dEdx2Tag"));
0139   }
0140 
0141   label_tv = consumes<TrackingVertexCollection>(pset.getParameter<edm::InputTag>("label_tv"));
0142   if (doPlotsOnlyForTruePV_ || doPVAssociationPlots_) {
0143     recoVertexToken_ = consumes<edm::View<reco::Vertex>>(pset.getUntrackedParameter<edm::InputTag>("label_vertex"));
0144     vertexAssociatorToken_ =
0145         consumes<reco::VertexToTrackingVertexAssociator>(pset.getUntrackedParameter<edm::InputTag>("vertexAssociator"));
0146   }
0147 
0148   if (doMVAPlots_) {
0149     mvaQualityCollectionTokens_.resize(labelToken.size());
0150     auto mvaPSet = pset.getUntrackedParameter<edm::ParameterSet>("mvaLabels");
0151     for (size_t iIter = 0; iIter < labelToken.size(); ++iIter) {
0152       edm::EDConsumerBase::Labels labels;
0153       labelsForToken(labelToken[iIter], labels);
0154       if (mvaPSet.exists(labels.module)) {
0155         mvaQualityCollectionTokens_[iIter] = edm::vector_transform(
0156             mvaPSet.getUntrackedParameter<std::vector<std::string>>(labels.module), [&](const std::string& tag) {
0157               return std::make_tuple(consumes<MVACollection>(edm::InputTag(tag, "MVAValues")),
0158                                      consumes<QualityMaskCollection>(edm::InputTag(tag, "QualityMasks")));
0159             });
0160       }
0161     }
0162   }
0163 
0164   tpSelector = TrackingParticleSelector(pset.getParameter<double>("ptMinTP"),
0165                                         pset.getParameter<double>("ptMaxTP"),
0166                                         pset.getParameter<double>("minRapidityTP"),
0167                                         pset.getParameter<double>("maxRapidityTP"),
0168                                         pset.getParameter<double>("tipTP"),
0169                                         pset.getParameter<double>("lipTP"),
0170                                         pset.getParameter<int>("minHitTP"),
0171                                         pset.getParameter<bool>("signalOnlyTP"),
0172                                         pset.getParameter<bool>("intimeOnlyTP"),
0173                                         pset.getParameter<bool>("chargedOnlyTP"),
0174                                         pset.getParameter<bool>("stableOnlyTP"),
0175                                         pset.getParameter<std::vector<int>>("pdgIdTP"),
0176                                         pset.getParameter<bool>("invertRapidityCutTP"),
0177                                         pset.getParameter<double>("minPhi"),
0178                                         pset.getParameter<double>("maxPhi"));
0179 
0180   cosmictpSelector = CosmicTrackingParticleSelector(pset.getParameter<double>("ptMinTP"),
0181                                                     pset.getParameter<double>("minRapidityTP"),
0182                                                     pset.getParameter<double>("maxRapidityTP"),
0183                                                     pset.getParameter<double>("tipTP"),
0184                                                     pset.getParameter<double>("lipTP"),
0185                                                     pset.getParameter<int>("minHitTP"),
0186                                                     pset.getParameter<bool>("chargedOnlyTP"),
0187                                                     pset.getParameter<std::vector<int>>("pdgIdTP"));
0188 
0189   ParameterSet psetVsPhi = psetForHistoProducerAlgo.getParameter<ParameterSet>("TpSelectorForEfficiencyVsPhi");
0190   dRtpSelector = TrackingParticleSelector(psetVsPhi.getParameter<double>("ptMin"),
0191                                           psetVsPhi.getParameter<double>("ptMax"),
0192                                           psetVsPhi.getParameter<double>("minRapidity"),
0193                                           psetVsPhi.getParameter<double>("maxRapidity"),
0194                                           psetVsPhi.getParameter<double>("tip"),
0195                                           psetVsPhi.getParameter<double>("lip"),
0196                                           psetVsPhi.getParameter<int>("minHit"),
0197                                           psetVsPhi.getParameter<bool>("signalOnly"),
0198                                           psetVsPhi.getParameter<bool>("intimeOnly"),
0199                                           psetVsPhi.getParameter<bool>("chargedOnly"),
0200                                           psetVsPhi.getParameter<bool>("stableOnly"),
0201                                           psetVsPhi.getParameter<std::vector<int>>("pdgId"),
0202                                           psetVsPhi.getParameter<bool>("invertRapidityCut"),
0203                                           psetVsPhi.getParameter<double>("minPhi"),
0204                                           psetVsPhi.getParameter<double>("maxPhi"));
0205 
0206   dRTrackSelector = MTVHistoProducerAlgoForTracker::makeRecoTrackSelectorFromTPSelectorParameters(psetVsPhi);
0207 
0208   useGsf = pset.getParameter<bool>("useGsf");
0209 
0210   _simHitTpMapTag = mayConsume<SimHitTPAssociationProducer::SimHitTPAssociationList>(
0211       pset.getParameter<edm::InputTag>("simHitTpMapTag"));
0212 
0213   if (calculateDrSingleCollection_) {
0214     labelTokenForDrCalculation =
0215         consumes<edm::View<reco::Track>>(pset.getParameter<edm::InputTag>("trackCollectionForDrCalculation"));
0216   }
0217 
0218   if (useAssociators_) {
0219     for (auto const& src : associators) {
0220       associatorTokens.push_back(consumes<reco::TrackToTrackingParticleAssociator>(src));
0221     }
0222   } else {
0223     for (auto const& src : associators) {
0224       associatormapStRs.push_back(consumes<reco::SimToRecoCollection>(src));
0225       associatormapRtSs.push_back(consumes<reco::RecoToSimCollection>(src));
0226     }
0227   }
0228 }
0229 
0230 MultiTrackValidator::~MultiTrackValidator() {}
0231 
0232 void MultiTrackValidator::bookHistograms(DQMStore::IBooker& ibook,
0233                                          edm::Run const&,
0234                                          edm::EventSetup const& setup,
0235                                          Histograms& histograms) const {
0236   if (label.empty()) {
0237     // Disable histogram booking if there are no track collections
0238     return;
0239   }
0240 
0241   const auto minColl = -0.5;
0242   const auto maxColl = label.size() - 0.5;
0243   const auto nintColl = label.size();
0244 
0245   auto binLabels = [&](dqm::reco::MonitorElement* me) {
0246     for (size_t i = 0; i < label.size(); ++i) {
0247       me->setBinLabel(i + 1, label[i].label());
0248     }
0249     me->disableAlphanumeric();
0250     return me;
0251   };
0252 
0253   //Booking histograms concerning with simulated tracks
0254   if (doSimPlots_) {
0255     ibook.cd();
0256     ibook.setCurrentFolder(dirName_ + "simulation");
0257 
0258     histoProducerAlgo_->bookSimHistos(ibook, histograms.histoProducerAlgo);
0259 
0260     ibook.cd();
0261     ibook.setCurrentFolder(dirName_);
0262   }
0263 
0264   for (unsigned int ww = 0; ww < associators.size(); ww++) {
0265     ibook.cd();
0266     // FIXME: these need to be moved to a subdirectory whose name depends on the associator
0267     ibook.setCurrentFolder(dirName_);
0268 
0269     if (doSummaryPlots_) {
0270       if (doSimTrackPlots_) {
0271         histograms.h_assoc_coll.push_back(
0272             binLabels(ibook.book1D("num_assoc(simToReco)_coll",
0273                                    "N of associated (simToReco) tracks vs track collection",
0274                                    nintColl,
0275                                    minColl,
0276                                    maxColl)));
0277         histograms.h_simul_coll.push_back(binLabels(
0278             ibook.book1D("num_simul_coll", "N of simulated tracks vs track collection", nintColl, minColl, maxColl)));
0279       }
0280       if (doRecoTrackPlots_) {
0281         histograms.h_reco_coll.push_back(binLabels(
0282             ibook.book1D("num_reco_coll", "N of reco track vs track collection", nintColl, minColl, maxColl)));
0283         histograms.h_assoc2_coll.push_back(
0284             binLabels(ibook.book1D("num_assoc(recoToSim)_coll",
0285                                    "N of associated (recoToSim) tracks vs track collection",
0286                                    nintColl,
0287                                    minColl,
0288                                    maxColl)));
0289         histograms.h_looper_coll.push_back(
0290             binLabels(ibook.book1D("num_duplicate_coll",
0291                                    "N of associated (recoToSim) looper tracks vs track collection",
0292                                    nintColl,
0293                                    minColl,
0294                                    maxColl)));
0295         histograms.h_pileup_coll.push_back(
0296             binLabels(ibook.book1D("num_pileup_coll",
0297                                    "N of associated (recoToSim) pileup tracks vs track collection",
0298                                    nintColl,
0299                                    minColl,
0300                                    maxColl)));
0301       }
0302     }
0303 
0304     for (unsigned int www = 0; www < label.size(); www++) {
0305       ibook.cd();
0306       InputTag algo = label[www];
0307       string dirName = dirName_;
0308       if (!algo.process().empty())
0309         dirName += algo.process() + "_";
0310       if (!algo.label().empty())
0311         dirName += algo.label() + "_";
0312       if (!algo.instance().empty())
0313         dirName += algo.instance() + "_";
0314       if (dirName.find("Tracks") < dirName.length()) {
0315         dirName.replace(dirName.find("Tracks"), 6, "");
0316       }
0317       string assoc = associators[ww].label();
0318       if (assoc.find("Track") < assoc.length()) {
0319         assoc.replace(assoc.find("Track"), 5, "");
0320       }
0321       dirName += assoc;
0322       std::replace(dirName.begin(), dirName.end(), ':', '_');
0323 
0324       ibook.setCurrentFolder(dirName);
0325 
0326       const bool doResolutionPlots = doResolutionPlots_[www];
0327 
0328       if (doSimTrackPlots_) {
0329         histoProducerAlgo_->bookSimTrackHistos(ibook, histograms.histoProducerAlgo, doResolutionPlots);
0330         if (doPVAssociationPlots_)
0331           histoProducerAlgo_->bookSimTrackPVAssociationHistos(ibook, histograms.histoProducerAlgo);
0332       }
0333 
0334       //Booking histograms concerning with reconstructed tracks
0335       if (doRecoTrackPlots_) {
0336         histoProducerAlgo_->bookRecoHistos(ibook, histograms.histoProducerAlgo, doResolutionPlots);
0337         if (dodEdxPlots_)
0338           histoProducerAlgo_->bookRecodEdxHistos(ibook, histograms.histoProducerAlgo);
0339         if (doPVAssociationPlots_)
0340           histoProducerAlgo_->bookRecoPVAssociationHistos(ibook, histograms.histoProducerAlgo);
0341         if (doMVAPlots_)
0342           histoProducerAlgo_->bookMVAHistos(
0343               ibook, histograms.histoProducerAlgo, mvaQualityCollectionTokens_[www].size());
0344       }
0345 
0346       if (doSeedPlots_) {
0347         histoProducerAlgo_->bookSeedHistos(ibook, histograms.histoProducerAlgo);
0348       }
0349     }  //end loop www
0350   }    // end loop ww
0351 }
0352 
0353 #ifdef EDM_ML_DEBUG
0354 namespace {
0355   void ensureEffIsSubsetOfFake(const TrackingParticleRefVector& eff, const TrackingParticleRefVector& fake) {
0356     // If efficiency RefVector is empty, don't check the product ids
0357     // as it will be 0:0 for empty. This covers also the case where
0358     // both are empty. The case of fake being empty and eff not is an
0359     // error.
0360     if (eff.empty())
0361       return;
0362 
0363     // First ensure product ids
0364     if (eff.id() != fake.id()) {
0365       throw cms::Exception("Configuration")
0366           << "Efficiency and fake TrackingParticle (refs) point to different collections (eff " << eff.id() << " fake "
0367           << fake.id()
0368           << "). This is not supported. Efficiency TP set must be the same or a subset of the fake TP set.";
0369     }
0370 
0371     // Same technique as in associationMapFilterValues
0372     edm::IndexSet fakeKeys;
0373     fakeKeys.reserve(fake.size());
0374     for (const auto& ref : fake) {
0375       fakeKeys.insert(ref.key());
0376     }
0377 
0378     for (const auto& ref : eff) {
0379       if (!fakeKeys.has(ref.key())) {
0380         throw cms::Exception("Configuration") << "Efficiency TrackingParticle " << ref.key()
0381                                               << " is not found from the set of fake TPs. This is not supported. The "
0382                                                  "efficiency TP set must be the same or a subset of the fake TP set.";
0383       }
0384     }
0385   }
0386 }  // namespace
0387 #endif
0388 
0389 const TrackingVertex::LorentzVector* MultiTrackValidator::getSimPVPosition(
0390     const edm::Handle<TrackingVertexCollection>& htv) const {
0391   for (const auto& simV : *htv) {
0392     if (simV.eventId().bunchCrossing() != 0)
0393       continue;  // remove OOTPU
0394     if (simV.eventId().event() != 0)
0395       continue;  // pick the PV of hard scatter
0396     return &(simV.position());
0397   }
0398   return nullptr;
0399 }
0400 
0401 const reco::Vertex::Point* MultiTrackValidator::getRecoPVPosition(
0402     const edm::Event& event, const edm::Handle<TrackingVertexCollection>& htv) const {
0403   edm::Handle<edm::View<reco::Vertex>> hvertex;
0404   event.getByToken(recoVertexToken_, hvertex);
0405 
0406   edm::Handle<reco::VertexToTrackingVertexAssociator> hvassociator;
0407   event.getByToken(vertexAssociatorToken_, hvassociator);
0408 
0409   auto v_r2s = hvassociator->associateRecoToSim(hvertex, htv);
0410   auto pvPtr = hvertex->refAt(0);
0411   if (pvPtr->isFake() || pvPtr->ndof() < 0)  // skip junk vertices
0412     return nullptr;
0413 
0414   auto pvFound = v_r2s.find(pvPtr);
0415   if (pvFound == v_r2s.end())
0416     return nullptr;
0417 
0418   for (const auto& vertexRefQuality : pvFound->val) {
0419     const TrackingVertex& tv = *(vertexRefQuality.first);
0420     if (tv.eventId().event() == 0 && tv.eventId().bunchCrossing() == 0) {
0421       return &(pvPtr->position());
0422     }
0423   }
0424 
0425   return nullptr;
0426 }
0427 
0428 void MultiTrackValidator::tpParametersAndSelection(
0429     const Histograms& histograms,
0430     const TrackingParticleRefVector& tPCeff,
0431     const edm::Event& event,
0432     const edm::EventSetup& setup,
0433     const reco::BeamSpot& bs,
0434     std::vector<std::tuple<TrackingParticle::Vector, TrackingParticle::Point>>& momVert_tPCeff,
0435     std::vector<size_t>& selected_tPCeff) const {
0436   selected_tPCeff.reserve(tPCeff.size());
0437   momVert_tPCeff.reserve(tPCeff.size());
0438   int nIntimeTPs = 0;
0439   if (parametersDefinerIsCosmic_) {
0440     for (size_t j = 0; j < tPCeff.size(); ++j) {
0441       const TrackingParticleRef& tpr = tPCeff[j];
0442 
0443       auto const& rec = parametersDefinerTP_->momentumAndVertex(event, setup, tpr);
0444       TrackingParticle::Vector const& momentum = std::get<0>(rec);
0445       TrackingParticle::Point const& vertex = std::get<1>(rec);
0446       if (doSimPlots_) {
0447         histoProducerAlgo_->fill_generic_simTrack_histos(
0448             histograms.histoProducerAlgo, momentum, vertex, tpr->eventId().bunchCrossing());
0449       }
0450       if (tpr->eventId().bunchCrossing() == 0)
0451         ++nIntimeTPs;
0452 
0453       if (cosmictpSelector(tpr, &bs, event, setup)) {
0454         selected_tPCeff.push_back(j);
0455         momVert_tPCeff.emplace_back(momentum, vertex);
0456       }
0457     }
0458   } else {
0459     size_t j = 0;
0460     for (auto const& tpr : tPCeff) {
0461       const TrackingParticle& tp = *tpr;
0462 
0463       // TODO: do we want to fill these from all TPs that include IT
0464       // and OOT (as below), or limit to IT+OOT TPs passing tpSelector
0465       // (as it was before)? The latter would require another instance
0466       // of tpSelector with intimeOnly=False.
0467       if (doSimPlots_) {
0468         histoProducerAlgo_->fill_generic_simTrack_histos(
0469             histograms.histoProducerAlgo, tp.momentum(), tp.vertex(), tp.eventId().bunchCrossing());
0470       }
0471       if (tp.eventId().bunchCrossing() == 0)
0472         ++nIntimeTPs;
0473 
0474       if (tpSelector(tp)) {
0475         selected_tPCeff.push_back(j);
0476         momVert_tPCeff.emplace_back(parametersDefinerTP_->momentumAndVertex(event, setup, tpr));
0477       }
0478       ++j;
0479     }
0480   }
0481   if (doSimPlots_) {
0482     histoProducerAlgo_->fill_simTrackBased_histos(histograms.histoProducerAlgo, nIntimeTPs);
0483   }
0484 }
0485 
0486 size_t MultiTrackValidator::tpDR(const TrackingParticleRefVector& tPCeff,
0487                                  const std::vector<size_t>& selected_tPCeff,
0488                                  DynArray<float>& dR_tPCeff,
0489                                  DynArray<float>& dR_tPCeff_jet,
0490                                  const edm::View<reco::Candidate>* cores) const {
0491   if (tPCeff.empty()) {
0492     return 0;
0493   }
0494   float etaL[tPCeff.size()], phiL[tPCeff.size()];
0495   size_t n_selTP_dr = 0;
0496   for (size_t iTP : selected_tPCeff) {
0497     //calculare dR wrt inclusive collection (also with PU, low pT, displaced)
0498     auto const& tp2 = *(tPCeff[iTP]);
0499     auto&& p = tp2.momentum();
0500     etaL[iTP] = etaFromXYZ(p.x(), p.y(), p.z());
0501     phiL[iTP] = atan2f(p.y(), p.x());
0502   }
0503   for (size_t iTP1 : selected_tPCeff) {
0504     auto const& tp = *(tPCeff[iTP1]);
0505     double dR = std::numeric_limits<double>::max();
0506     double dR_jet = std::numeric_limits<double>::max();
0507     if (dRtpSelector(tp)) {  //only for those needed for efficiency!
0508       ++n_selTP_dr;
0509       float eta = etaL[iTP1];
0510       float phi = phiL[iTP1];
0511       for (size_t iTP2 : selected_tPCeff) {
0512         //calculare dR wrt inclusive collection (also with PU, low pT, displaced)
0513         if (iTP1 == iTP2) {
0514           continue;
0515         }
0516         auto dR_tmp = reco::deltaR2(eta, phi, etaL[iTP2], phiL[iTP2]);
0517         if (dR_tmp < dR)
0518           dR = dR_tmp;
0519       }  // ttp2 (iTP)
0520       if (cores != nullptr) {
0521         for (unsigned int ji = 0; ji < cores->size(); ji++) {  //jet loop
0522           const reco::Candidate& jet = (*cores)[ji];
0523           double jet_eta = jet.eta();
0524           double jet_phi = jet.phi();
0525           auto dR_jet_tmp = reco::deltaR2(eta, phi, jet_eta, jet_phi);
0526           if (dR_jet_tmp < dR_jet)
0527             dR_jet = dR_jet_tmp;
0528         }
0529       }
0530     }
0531     dR_tPCeff[iTP1] = std::sqrt(dR);
0532     dR_tPCeff_jet[iTP1] = std::sqrt(dR_jet);
0533 
0534   }  // tp
0535   return n_selTP_dr;
0536 }
0537 
0538 void MultiTrackValidator::trackDR(const edm::View<reco::Track>& trackCollection,
0539                                   const edm::View<reco::Track>& trackCollectionDr,
0540                                   DynArray<float>& dR_trk,
0541                                   DynArray<float>& dR_trk_jet,
0542                                   const edm::View<reco::Candidate>* cores) const {
0543   if (trackCollectionDr.empty()) {
0544     return;
0545   }
0546   int i = 0;
0547   float etaL[trackCollectionDr.size()];
0548   float phiL[trackCollectionDr.size()];
0549   bool validL[trackCollectionDr.size()];
0550   for (auto const& track2 : trackCollectionDr) {
0551     auto&& p = track2.momentum();
0552     etaL[i] = etaFromXYZ(p.x(), p.y(), p.z());
0553     phiL[i] = atan2f(p.y(), p.x());
0554     validL[i] = !trackFromSeedFitFailed(track2);
0555     ++i;
0556   }
0557   for (View<reco::Track>::size_type i = 0; i < trackCollection.size(); ++i) {
0558     auto const& track = trackCollection[i];
0559     auto dR = std::numeric_limits<float>::max();
0560     auto dR_jet = std::numeric_limits<float>::max();
0561     if (!trackFromSeedFitFailed(track)) {
0562       auto&& p = track.momentum();
0563       float eta = etaFromXYZ(p.x(), p.y(), p.z());
0564       float phi = atan2f(p.y(), p.x());
0565       for (View<reco::Track>::size_type j = 0; j < trackCollectionDr.size(); ++j) {
0566         if (!validL[j])
0567           continue;
0568         auto dR_tmp = reco::deltaR2(eta, phi, etaL[j], phiL[j]);
0569         if ((dR_tmp < dR) & (dR_tmp > std::numeric_limits<float>::min()))
0570           dR = dR_tmp;
0571       }
0572       if (cores != nullptr) {
0573         for (unsigned int ji = 0; ji < cores->size(); ji++) {  //jet loop
0574           const reco::Candidate& jet = (*cores)[ji];
0575           double jet_eta = jet.eta();
0576           double jet_phi = jet.phi();
0577           auto dR_jet_tmp = reco::deltaR2(eta, phi, jet_eta, jet_phi);
0578           if (dR_jet_tmp < dR_jet)
0579             dR_jet = dR_jet_tmp;
0580         }
0581       }
0582     }
0583     dR_trk[i] = std::sqrt(dR);
0584     dR_trk_jet[i] = std::sqrt(dR_jet);
0585   }
0586 }
0587 
0588 void MultiTrackValidator::dqmAnalyze(const edm::Event& event,
0589                                      const edm::EventSetup& setup,
0590                                      const Histograms& histograms) const {
0591   if (label.empty()) {
0592     // Disable if there are no track collections
0593     return;
0594   }
0595 
0596   using namespace reco;
0597 
0598   LogDebug("TrackValidator") << "\n===================================================="
0599                              << "\n"
0600                              << "Analyzing new event"
0601                              << "\n"
0602                              << "====================================================\n"
0603                              << "\n";
0604 
0605   const TrackerTopology& ttopo = setup.getData(tTopoEsToken);
0606 
0607   // FIXME: we really need to move to edm::View for reading the
0608   // TrackingParticles... Unfortunately it has non-trivial
0609   // consequences on the associator/association interfaces etc.
0610   TrackingParticleRefVector tmpTPeff;
0611   TrackingParticleRefVector tmpTPfake;
0612   const TrackingParticleRefVector* tmpTPeffPtr = nullptr;
0613   const TrackingParticleRefVector* tmpTPfakePtr = nullptr;
0614 
0615   edm::Handle<TrackingParticleCollection> TPCollectionHeff;
0616   edm::Handle<TrackingParticleRefVector> TPCollectionHeffRefVector;
0617 
0618   const bool tp_effic_refvector = label_tp_effic.isUninitialized();
0619   if (!tp_effic_refvector) {
0620     event.getByToken(label_tp_effic, TPCollectionHeff);
0621     tmpTPeff.reserve(TPCollectionHeff->size());
0622     for (size_t i = 0, size = TPCollectionHeff->size(); i < size; ++i) {
0623       tmpTPeff.push_back(TrackingParticleRef(TPCollectionHeff, i));
0624     }
0625     tmpTPeffPtr = &tmpTPeff;
0626   } else {
0627     event.getByToken(label_tp_effic_refvector, TPCollectionHeffRefVector);
0628     tmpTPeffPtr = TPCollectionHeffRefVector.product();
0629   }
0630   if (!label_tp_fake.isUninitialized()) {
0631     edm::Handle<TrackingParticleCollection> TPCollectionHfake;
0632     event.getByToken(label_tp_fake, TPCollectionHfake);
0633     tmpTPfake.reserve(TPCollectionHfake->size());
0634     for (size_t i = 0, size = TPCollectionHfake->size(); i < size; ++i) {
0635       tmpTPfake.push_back(TrackingParticleRef(TPCollectionHfake, i));
0636     }
0637     tmpTPfakePtr = &tmpTPfake;
0638   } else {
0639     edm::Handle<TrackingParticleRefVector> TPCollectionHfakeRefVector;
0640     event.getByToken(label_tp_fake_refvector, TPCollectionHfakeRefVector);
0641     tmpTPfakePtr = TPCollectionHfakeRefVector.product();
0642   }
0643 
0644   TrackingParticleRefVector const& tPCeff = *tmpTPeffPtr;
0645   TrackingParticleRefVector const& tPCfake = *tmpTPfakePtr;
0646 
0647 #ifdef EDM_ML_DEBUG
0648   ensureEffIsSubsetOfFake(tPCeff, tPCfake);
0649 #endif
0650 
0651   if (parametersDefinerIsCosmic_) {
0652     edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;
0653     //warning: make sure the TP collection used in the map is the same used in the MTV!
0654     event.getByToken(_simHitTpMapTag, simHitsTPAssoc);
0655     parametersDefinerTP_->initEvent(simHitsTPAssoc);
0656     cosmictpSelector.initEvent(simHitsTPAssoc);
0657   }
0658 
0659   // Find the sim PV and tak its position
0660   edm::Handle<TrackingVertexCollection> htv;
0661   event.getByToken(label_tv, htv);
0662   const TrackingVertex::LorentzVector* theSimPVPosition = getSimPVPosition(htv);
0663   if (simPVMaxZ_ >= 0) {
0664     if (!theSimPVPosition)
0665       return;
0666     if (std::abs(theSimPVPosition->z()) > simPVMaxZ_)
0667       return;
0668   }
0669 
0670   // Check, when necessary, if reco PV matches to sim PV
0671   const reco::Vertex::Point* thePVposition = nullptr;
0672   if (doPlotsOnlyForTruePV_ || doPVAssociationPlots_) {
0673     thePVposition = getRecoPVPosition(event, htv);
0674     if (doPlotsOnlyForTruePV_ && !thePVposition)
0675       return;
0676 
0677     // Rest of the code assumes that if thePVposition is non-null, the
0678     // PV-association histograms get filled. In above, the "nullness"
0679     // is used to deliver the information if the reco PV is matched to
0680     // the sim PV.
0681     if (!doPVAssociationPlots_)
0682       thePVposition = nullptr;
0683   }
0684 
0685   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0686   event.getByToken(bsSrc, recoBeamSpotHandle);
0687   reco::BeamSpot const& bs = *recoBeamSpotHandle;
0688 
0689   edm::Handle<std::vector<PileupSummaryInfo>> puinfoH;
0690   event.getByToken(label_pileupinfo, puinfoH);
0691   PileupSummaryInfo puinfo;
0692 
0693   for (unsigned int puinfo_ite = 0; puinfo_ite < (*puinfoH).size(); ++puinfo_ite) {
0694     if ((*puinfoH)[puinfo_ite].getBunchCrossing() == 0) {
0695       puinfo = (*puinfoH)[puinfo_ite];
0696       break;
0697     }
0698   }
0699 
0700   // Number of 3D layers for TPs
0701   edm::Handle<edm::ValueMap<unsigned int>> tpNLayersH;
0702   event.getByToken(tpNLayersToken_, tpNLayersH);
0703   const auto& nLayers_tPCeff = *tpNLayersH;
0704 
0705   event.getByToken(tpNPixelLayersToken_, tpNLayersH);
0706   const auto& nPixelLayers_tPCeff = *tpNLayersH;
0707 
0708   event.getByToken(tpNStripStereoLayersToken_, tpNLayersH);
0709   const auto& nStripMonoAndStereoLayers_tPCeff = *tpNLayersH;
0710 
0711   // Precalculate TP selection (for efficiency), and momentum and vertex wrt PCA
0712   //
0713   // TODO: ParametersDefinerForTP ESProduct needs to be changed to
0714   // EDProduct because of consumes.
0715   //
0716   // In principle, we could just precalculate the momentum and vertex
0717   // wrt PCA for all TPs for once and put that to the event. To avoid
0718   // repetitive calculations those should be calculated only once for
0719   // each TP. That would imply that we should access TPs via Refs
0720   // (i.e. View) in here, since, in general, the eff and fake TP
0721   // collections can be different (and at least HI seems to use that
0722   // feature). This would further imply that the
0723   // RecoToSimCollection/SimToRecoCollection should be changed to use
0724   // View<TP> instead of vector<TP>, and migrate everything.
0725   //
0726   // Or we could take only one input TP collection, and do another
0727   // TP-selection to obtain the "fake" collection like we already do
0728   // for "efficiency" TPs.
0729   std::vector<size_t> selected_tPCeff;
0730   std::vector<std::tuple<TrackingParticle::Vector, TrackingParticle::Point>> momVert_tPCeff;
0731   tpParametersAndSelection(histograms, tPCeff, event, setup, bs, momVert_tPCeff, selected_tPCeff);
0732 
0733   //calculate dR for TPs
0734   declareDynArray(float, tPCeff.size(), dR_tPCeff);
0735 
0736   //calculate dR_jet for TPs
0737   const edm::View<reco::Candidate>* coresVector = nullptr;
0738   if (not cores_.isUninitialized()) {
0739     Handle<edm::View<reco::Candidate>> cores;
0740     event.getByToken(cores_, cores);
0741     if (cores.isValid()) {
0742       coresVector = cores.product();
0743     }
0744   }
0745   declareDynArray(float, tPCeff.size(), dR_tPCeff_jet);
0746 
0747   size_t n_selTP_dr = tpDR(tPCeff, selected_tPCeff, dR_tPCeff, dR_tPCeff_jet, coresVector);
0748 
0749   edm::Handle<View<Track>> trackCollectionForDrCalculation;
0750   if (calculateDrSingleCollection_) {
0751     event.getByToken(labelTokenForDrCalculation, trackCollectionForDrCalculation);
0752   }
0753 
0754   // dE/dx
0755   // at some point this could be generalized, with a vector of tags and a corresponding vector of Handles
0756   // I'm writing the interface such to take vectors of ValueMaps
0757   std::vector<const edm::ValueMap<reco::DeDxData>*> v_dEdx;
0758   if (dodEdxPlots_) {
0759     edm::Handle<edm::ValueMap<reco::DeDxData>> dEdx1Handle;
0760     edm::Handle<edm::ValueMap<reco::DeDxData>> dEdx2Handle;
0761     event.getByToken(m_dEdx1Tag, dEdx1Handle);
0762     event.getByToken(m_dEdx2Tag, dEdx2Handle);
0763     v_dEdx.push_back(dEdx1Handle.product());
0764     v_dEdx.push_back(dEdx2Handle.product());
0765   }
0766 
0767   std::vector<const MVACollection*> mvaCollections;
0768   std::vector<const QualityMaskCollection*> qualityMaskCollections;
0769   std::vector<float> mvaValues;
0770 
0771   int w = 0;  //counter counting the number of sets of histograms
0772   for (unsigned int ww = 0; ww < associators.size(); ww++) {
0773     // run value filtering of recoToSim map already here as it depends only on the association, not track collection
0774     reco::SimToRecoCollection const* simRecCollPFull = nullptr;
0775     reco::RecoToSimCollection const* recSimCollP = nullptr;
0776     reco::RecoToSimCollection recSimCollL;
0777     if (!useAssociators_) {
0778       Handle<reco::SimToRecoCollection> simtorecoCollectionH;
0779       event.getByToken(associatormapStRs[ww], simtorecoCollectionH);
0780       simRecCollPFull = simtorecoCollectionH.product();
0781 
0782       Handle<reco::RecoToSimCollection> recotosimCollectionH;
0783       event.getByToken(associatormapRtSs[ww], recotosimCollectionH);
0784       recSimCollP = recotosimCollectionH.product();
0785 
0786       // We need to filter the associations of the fake-TrackingParticle
0787       // collection only from RecoToSim collection, otherwise the
0788       // RecoToSim histograms get false entries
0789       recSimCollL = associationMapFilterValues(*recSimCollP, tPCfake);
0790       recSimCollP = &recSimCollL;
0791     }
0792 
0793     for (unsigned int www = 0; www < label.size();
0794          www++, w++) {  // need to increment w here, since there are many continues in the loop body
0795       //
0796       //get collections from the event
0797       //
0798       edm::Handle<View<Track>> trackCollectionHandle;
0799       if (!event.getByToken(labelToken[www], trackCollectionHandle) && ignoremissingtkcollection_)
0800         continue;
0801       const edm::View<Track>& trackCollection = *trackCollectionHandle;
0802 
0803       reco::SimToRecoCollection const* simRecCollP = nullptr;
0804       reco::SimToRecoCollection simRecCollL;
0805 
0806       //associate tracks
0807       LogTrace("TrackValidator") << "Analyzing " << label[www] << " with " << associators[ww] << "\n";
0808       if (useAssociators_) {
0809         edm::Handle<reco::TrackToTrackingParticleAssociator> theAssociator;
0810         event.getByToken(associatorTokens[ww], theAssociator);
0811 
0812         // The associator interfaces really need to be fixed...
0813         edm::RefToBaseVector<reco::Track> trackRefs;
0814         // trackRefs.vectorHolder()->reserve(trackCollection.size());  NOT a good idea
0815         for (edm::View<Track>::size_type i = 0; i < trackCollection.size(); ++i) {
0816           trackRefs.push_back(trackCollection.refAt(i));
0817         }
0818 
0819         LogTrace("TrackValidator") << "Calling associateRecoToSim method"
0820                                    << "\n";
0821         recSimCollL = theAssociator->associateRecoToSim(trackRefs, tPCfake);
0822         recSimCollP = &recSimCollL;
0823         LogTrace("TrackValidator") << "Calling associateSimToReco method"
0824                                    << "\n";
0825         // It is necessary to do the association wrt. fake TPs,
0826         // because this SimToReco association is used also for
0827         // duplicates. Since the set of efficiency TPs are required to
0828         // be a subset of the set of fake TPs, for efficiency
0829         // histograms it doesn't matter if the association contains
0830         // associations of TPs not in the set of efficiency TPs.
0831         simRecCollL = theAssociator->associateSimToReco(trackRefs, tPCfake);
0832         simRecCollP = &simRecCollL;
0833       } else {
0834         // We need to filter the associations of the current track
0835         // collection only from SimToReco collection, otherwise the
0836         // SimToReco histograms get false entries. The filtering must
0837         // be done separately for each track collection.
0838         simRecCollL = associationMapFilterValues(*simRecCollPFull, trackCollection);
0839         simRecCollP = &simRecCollL;
0840       }
0841 
0842       reco::RecoToSimCollection const& recSimColl = *recSimCollP;
0843       reco::SimToRecoCollection const& simRecColl = *simRecCollP;
0844 
0845       // read MVA collections
0846       if (doMVAPlots_ && !mvaQualityCollectionTokens_[www].empty()) {
0847         edm::Handle<MVACollection> hmva;
0848         edm::Handle<QualityMaskCollection> hqual;
0849         for (const auto& tokenTpl : mvaQualityCollectionTokens_[www]) {
0850           event.getByToken(std::get<0>(tokenTpl), hmva);
0851           event.getByToken(std::get<1>(tokenTpl), hqual);
0852 
0853           mvaCollections.push_back(hmva.product());
0854           qualityMaskCollections.push_back(hqual.product());
0855           if (mvaCollections.back()->size() != trackCollection.size()) {
0856             throw cms::Exception("Configuration")
0857                 << "Inconsistency in track collection and MVA sizes. Track collection " << www << " has "
0858                 << trackCollection.size() << " tracks, whereas the MVA " << (mvaCollections.size() - 1)
0859                 << " for it has " << mvaCollections.back()->size() << " entries. Double-check your configuration.";
0860           }
0861           if (qualityMaskCollections.back()->size() != trackCollection.size()) {
0862             throw cms::Exception("Configuration")
0863                 << "Inconsistency in track collection and quality mask sizes. Track collection " << www << " has "
0864                 << trackCollection.size() << " tracks, whereas the quality mask " << (qualityMaskCollections.size() - 1)
0865                 << " for it has " << qualityMaskCollections.back()->size()
0866                 << " entries. Double-check your configuration.";
0867           }
0868         }
0869       }
0870 
0871       // ########################################################
0872       // fill simulation histograms (LOOP OVER TRACKINGPARTICLES)
0873       // ########################################################
0874 
0875       //compute number of tracks per eta interval
0876       //
0877       LogTrace("TrackValidator") << "\n# of TrackingParticles: " << tPCeff.size() << "\n";
0878       int st(0);  //This counter counts the number of simulated tracks passing the MTV selection (i.e. tpSelector(tp) )
0879 
0880       //loop over already-selected TPs for tracking efficiency
0881       for (size_t i = 0; i < selected_tPCeff.size(); ++i) {
0882         size_t iTP = selected_tPCeff[i];
0883         const TrackingParticleRef& tpr = tPCeff[iTP];
0884         const TrackingParticle& tp = *tpr;
0885 
0886         auto const& momVert = momVert_tPCeff[i];
0887         TrackingParticle::Vector momentumTP;
0888         TrackingParticle::Point vertexTP;
0889 
0890         double dxySim(0);
0891         double dzSim(0);
0892         double dxyPVSim = 0;
0893         double dzPVSim = 0;
0894         double dR = dR_tPCeff[iTP];
0895         double dR_jet = dR_tPCeff_jet[iTP];
0896 
0897         //---------- THIS PART HAS TO BE CLEANED UP. THE PARAMETER DEFINER WAS NOT MEANT TO BE USED IN THIS WAY ----------
0898         //If the TrackingParticle is collison like, get the momentum and vertex at production state
0899         if (!parametersDefinerIsCosmic_) {
0900           momentumTP = tp.momentum();
0901           vertexTP = tp.vertex();
0902           //Calcualte the impact parameters w.r.t. PCA
0903           const TrackingParticle::Vector& momentum = std::get<TrackingParticle::Vector>(momVert);
0904           const TrackingParticle::Point& vertex = std::get<TrackingParticle::Point>(momVert);
0905           dxySim = TrackingParticleIP::dxy(vertex, momentum, bs.position());
0906           dzSim = TrackingParticleIP::dz(vertex, momentum, bs.position());
0907 
0908           if (theSimPVPosition) {
0909             dxyPVSim = TrackingParticleIP::dxy(vertex, momentum, *theSimPVPosition);
0910             dzPVSim = TrackingParticleIP::dz(vertex, momentum, *theSimPVPosition);
0911           }
0912         }
0913         //If the TrackingParticle is comics, get the momentum and vertex at PCA
0914         else {
0915           momentumTP = std::get<TrackingParticle::Vector>(momVert);
0916           vertexTP = std::get<TrackingParticle::Point>(momVert);
0917           dxySim = TrackingParticleIP::dxy(vertexTP, momentumTP, bs.position());
0918           dzSim = TrackingParticleIP::dz(vertexTP, momentumTP, bs.position());
0919 
0920           // Do dxy and dz vs. PV make any sense for cosmics? I guess not
0921         }
0922         //---------- THE PART ABOVE HAS TO BE CLEANED UP. THE PARAMETER DEFINER WAS NOT MEANT TO BE USED IN THIS WAY ----------
0923 
0924         // in the coming lines, histos are filled using as input
0925         // - momentumTP
0926         // - vertexTP
0927         // - dxySim
0928         // - dzSim
0929         if (!doSimTrackPlots_)
0930           continue;
0931 
0932         // ##############################################
0933         // fill RecoAssociated SimTracks' histograms
0934         // ##############################################
0935         const reco::Track* matchedTrackPointer = nullptr;
0936         const reco::Track* matchedSecondTrackPointer = nullptr;
0937         unsigned int selectsLoose = mvaCollections.size();
0938         unsigned int selectsHP = mvaCollections.size();
0939         if (simRecColl.find(tpr) != simRecColl.end()) {
0940           auto const& rt = simRecColl[tpr];
0941           if (!rt.empty()) {
0942             // isRecoMatched = true; // UNUSED
0943             matchedTrackPointer = rt.begin()->first.get();
0944             if (rt.size() >= 2) {
0945               matchedSecondTrackPointer = (rt.begin() + 1)->first.get();
0946             }
0947             LogTrace("TrackValidator") << "TrackingParticle #" << st << " with pt=" << sqrt(momentumTP.perp2())
0948                                        << " associated with quality:" << rt.begin()->second << "\n";
0949 
0950             if (doMVAPlots_) {
0951               // for each MVA we need to take the value of the track
0952               // with largest MVA value (for the cumulative histograms)
0953               //
0954               // also identify the first MVA that possibly selects any
0955               // track matched to this TrackingParticle, separately
0956               // for loose and highPurity qualities
0957               for (size_t imva = 0; imva < mvaCollections.size(); ++imva) {
0958                 const auto& mva = *(mvaCollections[imva]);
0959                 const auto& qual = *(qualityMaskCollections[imva]);
0960 
0961                 auto iMatch = rt.begin();
0962                 float maxMva = mva[iMatch->first.key()];
0963                 for (; iMatch != rt.end(); ++iMatch) {
0964                   auto itrk = iMatch->first.key();
0965                   maxMva = std::max(maxMva, mva[itrk]);
0966 
0967                   if (selectsLoose >= imva && trackSelected(qual[itrk], reco::TrackBase::loose))
0968                     selectsLoose = imva;
0969                   if (selectsHP >= imva && trackSelected(qual[itrk], reco::TrackBase::highPurity))
0970                     selectsHP = imva;
0971                 }
0972                 mvaValues.push_back(maxMva);
0973               }
0974             }
0975           }
0976         } else {
0977           LogTrace("TrackValidator") << "TrackingParticle #" << st << " with pt,eta,phi: " << sqrt(momentumTP.perp2())
0978                                      << " , " << momentumTP.eta() << " , " << momentumTP.phi() << " , "
0979                                      << " NOT associated to any reco::Track"
0980                                      << "\n";
0981         }
0982 
0983         int nSimHits = tp.numberOfTrackerHits();
0984         int nSimLayers = nLayers_tPCeff[tpr];
0985         int nSimPixelLayers = nPixelLayers_tPCeff[tpr];
0986         int nSimStripMonoAndStereoLayers = nStripMonoAndStereoLayers_tPCeff[tpr];
0987         histoProducerAlgo_->fill_recoAssociated_simTrack_histos(histograms.histoProducerAlgo,
0988                                                                 w,
0989                                                                 tp,
0990                                                                 momentumTP,
0991                                                                 vertexTP,
0992                                                                 dxySim,
0993                                                                 dzSim,
0994                                                                 dxyPVSim,
0995                                                                 dzPVSim,
0996                                                                 nSimHits,
0997                                                                 nSimLayers,
0998                                                                 nSimPixelLayers,
0999                                                                 nSimStripMonoAndStereoLayers,
1000                                                                 matchedTrackPointer,
1001                                                                 puinfo.getPU_NumInteractions(),
1002                                                                 dR,
1003                                                                 dR_jet,
1004                                                                 thePVposition,
1005                                                                 theSimPVPosition,
1006                                                                 bs.position(),
1007                                                                 mvaValues,
1008                                                                 selectsLoose,
1009                                                                 selectsHP);
1010         mvaValues.clear();
1011 
1012         if (matchedTrackPointer && matchedSecondTrackPointer) {
1013           histoProducerAlgo_->fill_duplicate_histos(
1014               histograms.histoProducerAlgo, w, *matchedTrackPointer, *matchedSecondTrackPointer);
1015         }
1016 
1017         if (doSummaryPlots_) {
1018           if (dRtpSelector(tp)) {
1019             histograms.h_simul_coll[ww]->Fill(www);
1020             if (matchedTrackPointer) {
1021               histograms.h_assoc_coll[ww]->Fill(www);
1022             }
1023           }
1024         }
1025 
1026       }  // End  for (TrackingParticleCollection::size_type i=0; i<tPCeff.size(); i++){
1027 
1028       // ##############################################
1029       // fill recoTracks histograms (LOOP OVER TRACKS)
1030       // ##############################################
1031       if (!doRecoTrackPlots_)
1032         continue;
1033       LogTrace("TrackValidator") << "\n# of reco::Tracks with " << label[www].process() << ":" << label[www].label()
1034                                  << ":" << label[www].instance() << ": " << trackCollection.size() << "\n";
1035 
1036       int at(0);  //This counter counts the number of recoTracks that are associated to SimTracks
1037       int rT(0);  //This counter counts the number of recoTracks in general
1038       int seed_fit_failed = 0;
1039       size_t n_selTrack_dr = 0;
1040 
1041       //calculate dR for tracks
1042       declareDynArray(float, trackCollection.size(), dR_trk);
1043       declareDynArray(float, trackCollection.size(), dR_trk_jet);
1044 #ifndef NO_TRACK_DR
1045       // this accounts for most of the time spent in MTV and it is used to fill just one histo that is of doubtful usefulness but (maybe) for the whole collection
1046       const edm::View<Track>* trackCollectionDr = &trackCollection;
1047       if (calculateDrSingleCollection_) {
1048         trackCollectionDr = trackCollectionForDrCalculation.product();
1049       }
1050       trackDR(trackCollection, *trackCollectionDr, dR_trk, dR_trk_jet, coresVector);
1051 #endif
1052 
1053       for (View<Track>::size_type i = 0; i < trackCollection.size(); ++i) {
1054         auto track = trackCollection.refAt(i);
1055         rT++;
1056         if (trackFromSeedFitFailed(*track))
1057           ++seed_fit_failed;
1058         if ((*dRTrackSelector)(*track, bs.position()))
1059           ++n_selTrack_dr;
1060 
1061         bool isSigSimMatched(false);
1062         bool isSimMatched(false);
1063         bool isChargeMatched(true);
1064         int numAssocRecoTracks = 0;
1065         int nSimHits = 0;
1066         double sharedFraction = 0.;
1067 
1068         auto tpFound = recSimColl.find(track);
1069         isSimMatched = tpFound != recSimColl.end();
1070         if (applyTPSelToSimMatch_ && isSimMatched)
1071           isSimMatched = tpSelector(*tpFound->val[0].first);
1072         if (isSimMatched) {
1073           const auto& tp = tpFound->val;
1074           nSimHits = tp[0].first->numberOfTrackerHits();
1075           sharedFraction = tp[0].second;
1076           if (tp[0].first->charge() != track->charge())
1077             isChargeMatched = false;
1078           if (simRecColl.find(tp[0].first) != simRecColl.end())
1079             numAssocRecoTracks = simRecColl[tp[0].first].size();
1080           at++;
1081           for (unsigned int tp_ite = 0; tp_ite < tp.size(); ++tp_ite) {
1082             TrackingParticle trackpart = *(tp[tp_ite].first);
1083             if ((trackpart.eventId().event() == 0) && (trackpart.eventId().bunchCrossing() == 0)) {
1084               isSigSimMatched = true;
1085               break;
1086             }
1087           }
1088           LogTrace("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
1089                                      << " associated with quality:" << tp.begin()->second << "\n";
1090         } else {
1091           LogTrace("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
1092                                      << " NOT associated to any TrackingParticle"
1093                                      << "\n";
1094         }
1095 
1096         // set MVA values for this track
1097         // take also the indices of first MVAs to select by loose and
1098         // HP quality
1099         unsigned int selectsLoose = mvaCollections.size();
1100         unsigned int selectsHP = mvaCollections.size();
1101         if (doMVAPlots_) {
1102           for (size_t imva = 0; imva < mvaCollections.size(); ++imva) {
1103             const auto& mva = *(mvaCollections[imva]);
1104             const auto& qual = *(qualityMaskCollections[imva]);
1105             mvaValues.push_back(mva[i]);
1106 
1107             if (selectsLoose >= imva && trackSelected(qual[i], reco::TrackBase::loose))
1108               selectsLoose = imva;
1109             if (selectsHP >= imva && trackSelected(qual[i], reco::TrackBase::highPurity))
1110               selectsHP = imva;
1111           }
1112         }
1113 
1114         double dR = dR_trk[i];
1115         double dR_jet = dR_trk_jet[i];
1116         histoProducerAlgo_->fill_generic_recoTrack_histos(histograms.histoProducerAlgo,
1117                                                           w,
1118                                                           *track,
1119                                                           ttopo,
1120                                                           bs.position(),
1121                                                           thePVposition,
1122                                                           theSimPVPosition,
1123                                                           isSimMatched,
1124                                                           isSigSimMatched,
1125                                                           isChargeMatched,
1126                                                           numAssocRecoTracks,
1127                                                           puinfo.getPU_NumInteractions(),
1128                                                           nSimHits,
1129                                                           sharedFraction,
1130                                                           dR,
1131                                                           dR_jet,
1132                                                           mvaValues,
1133                                                           selectsLoose,
1134                                                           selectsHP);
1135         mvaValues.clear();
1136 
1137         if (doSummaryPlots_) {
1138           histograms.h_reco_coll[ww]->Fill(www);
1139           if (isSimMatched) {
1140             histograms.h_assoc2_coll[ww]->Fill(www);
1141             if (numAssocRecoTracks > 1) {
1142               histograms.h_looper_coll[ww]->Fill(www);
1143             }
1144             if (!isSigSimMatched) {
1145               histograms.h_pileup_coll[ww]->Fill(www);
1146             }
1147           }
1148         }
1149 
1150         // dE/dx
1151         if (dodEdxPlots_)
1152           histoProducerAlgo_->fill_dedx_recoTrack_histos(histograms.histoProducerAlgo, w, track, v_dEdx);
1153 
1154         //Fill other histos
1155         if (!isSimMatched)
1156           continue;
1157 
1158         histoProducerAlgo_->fill_simAssociated_recoTrack_histos(histograms.histoProducerAlgo, w, *track);
1159 
1160         /* TO BE FIXED LATER
1161     if (associators[ww]=="trackAssociatorByChi2"){
1162       //association chi2
1163       double assocChi2 = -tp.begin()->second;//in association map is stored -chi2
1164       h_assochi2[www]->Fill(assocChi2);
1165       h_assochi2_prob[www]->Fill(TMath::Prob((assocChi2)*5,5));
1166     }
1167     else if (associators[ww]=="quickTrackAssociatorByHits"){
1168       double fraction = tp.begin()->second;
1169       h_assocFraction[www]->Fill(fraction);
1170       h_assocSharedHit[www]->Fill(fraction*track->numberOfValidHits());
1171     }
1172     */
1173 
1174         if (doResolutionPlots_[www]) {
1175           //Get tracking particle parameters at point of closest approach to the beamline
1176           TrackingParticleRef tpr = tpFound->val.begin()->first;
1177           TrackingParticle::Vector momentumTP = parametersDefinerTP_->momentum(event, setup, tpr);
1178           TrackingParticle::Point vertexTP = parametersDefinerTP_->vertex(event, setup, tpr);
1179           int chargeTP = tpr->charge();
1180 
1181           histoProducerAlgo_->fill_ResoAndPull_recoTrack_histos(
1182               histograms.histoProducerAlgo, w, momentumTP, vertexTP, chargeTP, *track, bs.position());
1183         }
1184 
1185         //TO BE FIXED
1186         //std::vector<PSimHit> simhits=tpr.get()->trackPSimHit(DetId::Tracker);
1187         //nrecHit_vs_nsimHit_rec2sim[w]->Fill(track->numberOfValidHits(), (int)(simhits.end()-simhits.begin() ));
1188 
1189       }  // End of for(View<Track>::size_type i=0; i<trackCollection.size(); ++i){
1190       mvaCollections.clear();
1191       qualityMaskCollections.clear();
1192 
1193       histoProducerAlgo_->fill_trackBased_histos(histograms.histoProducerAlgo, w, at, rT, n_selTrack_dr, n_selTP_dr);
1194       // Fill seed-specific histograms
1195       if (doSeedPlots_) {
1196         histoProducerAlgo_->fill_seed_histos(
1197             histograms.histoProducerAlgo, www, seed_fit_failed, trackCollection.size());
1198       }
1199 
1200       LogTrace("TrackValidator") << "Collection " << www << "\n"
1201                                  << "Total Simulated (selected): " << n_selTP_dr << "\n"
1202                                  << "Total Reconstructed (selected): " << n_selTrack_dr << "\n"
1203                                  << "Total Reconstructed: " << rT << "\n"
1204                                  << "Total Associated (recoToSim): " << at << "\n"
1205                                  << "Total Fakes: " << rT - at << "\n";
1206     }  // End of  for (unsigned int www=0;www<label.size();www++){
1207   }    //END of for (unsigned int ww=0;ww<associators.size();ww++){
1208 }
[ 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!