Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-09-16 03:24:22

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