VertexAssociation/src/VertexAssociatorByTracks.cc

0001 #include "SimTracker/VertexAssociation/interface/VertexAssociatorByTracks.h"
0002
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "DataFormats/Common/interface/Ref.h"
0005 #include "DataFormats/Common/interface/RefToBase.h"
0006 #include "DataFormats/VertexReco/interface/Vertex.h"
0007 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0008
0009 #include "FWCore/Framework/interface/Frameworkfwd.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
0013 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0014 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0015
0016 #include "SimTracker/Common/interface/TrackingParticleSelector.h"
0017
0018 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0019
0020 /* Constructor */
0021 VertexAssociatorByTracks::VertexAssociatorByTracks(const edm::EDProductGetter *productGetter,
0022                                                    double R2SMatchedSimRatio,
0023                                                    double R2SMatchedRecoRatio,
0024                                                    double S2RMatchedSimRatio,
0025                                                    double S2RMatchedRecoRatio,
0026                                                    const TrackingParticleSelector *selector,
0027                                                    reco::TrackBase::TrackQuality trackQuality,
0028                                                    const reco::RecoToSimCollection *trackRecoToSimAssociation,
0029                                                    const reco::SimToRecoCollection *trackSimToRecoAssociation)
0030     : productGetter_(productGetter),
0031       R2SMatchedSimRatio_(R2SMatchedSimRatio),
0032       R2SMatchedRecoRatio_(R2SMatchedRecoRatio),
0033       S2RMatchedSimRatio_(S2RMatchedSimRatio),
0034       S2RMatchedRecoRatio_(S2RMatchedRecoRatio),
0035       selector_(selector),
0036       trackQuality_(trackQuality),
0037       trackRecoToSimAssociation_(trackRecoToSimAssociation),
0038       trackSimToRecoAssociation_(trackSimToRecoAssociation) {}
0039
0040 /* Destructor */
0041 VertexAssociatorByTracks::~VertexAssociatorByTracks() {}
0042
0043 reco::VertexRecoToSimCollection VertexAssociatorByTracks::associateRecoToSim(
0044     const edm::Handle<edm::View<reco::Vertex>> &recoVertexes,
0045     const edm::Handle<TrackingVertexCollection> &trackingVertexes) const {
0046   reco::VertexRecoToSimCollection outputCollection(productGetter_);
0047
0048   std::map<TrackingVertexRef, std::pair<double, std::size_t>> matches;
0049
0050   LogDebug("VertexAssociation") << "reco::VertexCollection size = " << recoVertexes->size()
0051                                 << " ; TrackingVertexCollection size = " << trackingVertexes->size() << std::endl;
0052
0053   // Loop over RecoVertex
0054   for (std::size_t recoIndex = 0; recoIndex < recoVertexes->size(); ++recoIndex) {
0055     reco::VertexBaseRef recoVertex(recoVertexes, recoIndex);
0056
0057     matches.clear();
0058
0059     double recoDaughterWeight = 0.;
0060
0061     // Loop over daughter tracks of RecoVertex
0062     for (reco::Vertex::trackRef_iterator recoDaughter = recoVertex->tracks_begin();
0063          recoDaughter != recoVertex->tracks_end();
0064          ++recoDaughter) {
0065       // Check the quality of the RecoDoughters
0066       if (!(*recoDaughter)->quality(trackQuality_))
0067         continue;
0068
0069       // Check for association for the given RecoDaughter
0070       if (trackRecoToSimAssociation_->numberOfAssociations(*recoDaughter) > 0) {
0071         std::vector<std::pair<TrackingParticleRef, double>> associations = (*trackRecoToSimAssociation_)[*recoDaughter];
0072
0073         // Loop over TrackingParticles associated with RecoDaughter
0074         for (std::vector<std::pair<TrackingParticleRef, double>>::const_iterator association = associations.begin();
0075              association != associations.end();
0076              ++association) {
0077           // Get a reference to parent vertex of TrackingParticle associated to
0078           // the RecoDaughter
0079           TrackingVertexRef trackingVertex = association->first->parentVertex();
0080           // Store matched RecoDaughter to the trackingVertex
0081           matches[trackingVertex].first += recoVertex->trackWeight(*recoDaughter);
0082           matches[trackingVertex].second++;
0083         }
0084       }
0085
0086       recoDaughterWeight += recoVertex->trackWeight(*recoDaughter);
0087     }  // loop over daughter tracks of RecoVertex
0088
0089     std::size_t assoIndex = 0;
0090
0091     // Loop over map between TrackingVertexes and matched RecoDaugther
0092     for (std::map<TrackingVertexRef, std::pair<double, std::size_t>>::const_iterator match = matches.begin();
0093          match != matches.end();
0094          ++match) {
0095       // Getting the TrackingVertex information
0096       TrackingVertexRef trackingVertex = match->first;
0097       double matchedDaughterWeight = match->second.first;
0098       std::size_t matchedDaughterCounter = match->second.second;
0099
0100       // Count for only reconstructible SimDaughters
0101       std::size_t simDaughterCounter = 0;
0102
0103       for (TrackingVertex::tp_iterator simDaughter = trackingVertex->daughterTracks_begin();
0104            simDaughter != trackingVertex->daughterTracks_end();
0105            ++simDaughter)
0106         if ((*selector_)(**simDaughter))
0107           simDaughterCounter++;
0108
0109       // Sanity condition in case that reconstructable condition is too tight
0110       if (simDaughterCounter < matchedDaughterCounter)
0111         simDaughterCounter = matchedDaughterCounter;
0112
0113       // Condition over S2RMatchedSimRatio
0114       if ((double)matchedDaughterCounter / simDaughterCounter < R2SMatchedSimRatio_)
0115         continue;
0116
0117       double quality = (double)matchedDaughterWeight / recoDaughterWeight;
0118
0119       // Condition over R2SMatchedRecoRatio
0120       if (quality < R2SMatchedRecoRatio_)
0121         continue;
0122
0123       outputCollection.insert(recoVertex, std::make_pair(trackingVertex, quality));
0124
0125       LogTrace("VertexAssociation") << "R2S: INDEX " << assoIndex << " ; SimDaughterCounter = " << simDaughterCounter
0126                                     << " ; RecoDaughterWeight = " << recoDaughterWeight
0127                                     << " ; MatchedDaughterCounter = " << matchedDaughterCounter
0128                                     << " ; MatchedDaughterWeight = " << matchedDaughterWeight
0129                                     << " ; quality = " << quality;
0130
0131       assoIndex++;
0132     }
0133   }  // Loop on RecoVertex
0134
0135   LogTrace("VertexAssociation") << "\nRecoToSim OUTPUT COLLECTION: outputCollection.size() = "
0136                                 << outputCollection.size() << std::endl;
0137
0138   return outputCollection;
0139 }
0140
0141 reco::VertexSimToRecoCollection VertexAssociatorByTracks::associateSimToReco(
0142     const edm::Handle<edm::View<reco::Vertex>> &recoVertexes,
0143     const edm::Handle<TrackingVertexCollection> &trackingVertexes) const {
0144   reco::VertexSimToRecoCollection outputCollection(productGetter_);
0145
0146   // Loop over TrackingVertexes
0147   std::map<std::size_t, std::pair<double, std::size_t>> matches;
0148
0149   // Loop over TrackingVertexes
0150   for (std::size_t simIndex = 0; simIndex < trackingVertexes->size(); ++simIndex) {
0151     TrackingVertexRef trackingVertex(trackingVertexes, simIndex);
0152
0153     matches.clear();
0154
0155     std::size_t simDaughterCounter = 0;
0156
0157     // Loop over daughter tracks of TrackingVertex
0158     for (TrackingVertex::tp_iterator simDaughter = trackingVertex->daughterTracks_begin();
0159          simDaughter != trackingVertex->daughterTracks_end();
0160          ++simDaughter) {
0161       // Select only reconstructible SimDaughters
0162       if (!(*selector_)(**simDaughter))
0163         continue;
0164
0165       // Check for association for the given RecoDaughter
0166       if (trackSimToRecoAssociation_->numberOfAssociations(*simDaughter) > 0) {
0167         std::vector<std::pair<reco::TrackBaseRef, double>> associations = (*trackSimToRecoAssociation_)[*simDaughter];
0168
0169         // Loop over RecoTracks associated with TrackingParticle
0170         for (std::vector<std::pair<reco::TrackBaseRef, double>>::const_iterator association = associations.begin();
0171              association != associations.end();
0172              ++association) {
0173           reco::TrackBaseRef recoTrack = association->first;
0174
0175           for (std::size_t recoIndex = 0; recoIndex < recoVertexes->size(); ++recoIndex) {
0176             reco::VertexBaseRef recoVertex(recoVertexes, recoIndex);
0177
0178             for (reco::Vertex::trackRef_iterator recoDaughter = recoVertex->tracks_begin();
0179                  recoDaughter != recoVertex->tracks_end();
0180                  ++recoDaughter) {
0181               // Store matched RecoDaughter to the RecoVertex
0182               if (recoDaughter->id() == recoTrack.id() && recoDaughter->key() == recoTrack.key()) {
0183                 matches[recoIndex].first += recoVertex->trackWeight(*recoDaughter);
0184                 matches[recoIndex].second++;
0185               }
0186             }
0187           }
0188         }  // loop a recotracks
0189       }
0190
0191       simDaughterCounter++;
0192     }  // loop a simDaughter
0193
0194     std::size_t assoIndex = 0;
0195
0196     // Loop over map, set score, add to outputCollection
0197     for (std::map<std::size_t, std::pair<double, std::size_t>>::const_iterator match = matches.begin();
0198          match != matches.end();
0199          ++match) {
0200       // Getting the TrackingVertex information
0201       reco::VertexBaseRef recoVertex(recoVertexes, match->first);
0202       double matchedDaughterWeight = match->second.first;
0203       std::size_t matchedDaughterCounter = match->second.second;
0204
0205       double recoDaughterWeight = 0.;
0206
0207       // Weighted count of those tracks with a given quality of the
0208       // RecoDoughters
0209       for (reco::Vertex::trackRef_iterator recoDaughter = recoVertex->tracks_begin();
0210            recoDaughter != recoVertex->tracks_end();
0211            ++recoDaughter)
0212         if ((*recoDaughter)->quality(trackQuality_))
0213           recoDaughterWeight += recoVertex->trackWeight(*recoDaughter);
0214
0215       // Sanity condition in case that track quality condition is too tight
0216       if (recoDaughterWeight < matchedDaughterWeight)
0217         recoDaughterWeight = matchedDaughterWeight;
0218
0219       // Condition over S2RMatchedRecoRatio
0220       if (matchedDaughterWeight / recoDaughterWeight < S2RMatchedRecoRatio_)
0221         continue;
0222
0223       double quality = (double)matchedDaughterCounter / simDaughterCounter;
0224
0225       // Condition over S2RMatchedSimRatio
0226       if (quality < S2RMatchedSimRatio_)
0227         continue;
0228
0229       outputCollection.insert(trackingVertex, std::make_pair(recoVertex, quality));
0230
0231       LogTrace("VertexAssociation") << "R2S: INDEX " << assoIndex << " ; SimDaughterCounter = " << simDaughterCounter
0232                                     << " ; RecoDaughterWeight = " << recoDaughterWeight
0233                                     << " ; MatchedDaughterCounter = " << matchedDaughterCounter
0234                                     << " ; MatchedDaughterWeight = " << matchedDaughterWeight
0235                                     << " ; quality = " << quality;
0236
0237       assoIndex++;
0238     }
0239   }  // loop over TrackingVertexes
0240
0241   LogTrace("VertexAssociation") << "SimToReco OUTPUT COLLECTION: outputCollection.size() = " << outputCollection.size()
0242                                 << std::endl;
0243
0244   return outputCollection;
0245 }