RecoTrack/plugins/SimpleTrackValidation.cc

0001 // system include files
0002 #include <memory>
0003
0004 #include "TTree.h"
0005 #include "TFile.h"
0006
0007 // user include files
0008 #include "FWCore/Framework/interface/Frameworkfwd.h"
0009 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0010
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/Utilities/interface/InputTag.h"
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0018
0019 #include "SimTracker/TrackerHitAssociation/interface/ClusterTPAssociation.h"
0020 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0021 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0022
0023 #include "SimTracker/Common/interface/TrackingParticleSelector.h"
0024
0025 #include "FWCore/ServiceRegistry/interface/Service.h"
0026 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0027
0028 using reco::TrackCollection;
0029
0030 class SimpleTrackValidation : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0031 public:
0032   explicit SimpleTrackValidation(const edm::ParameterSet&);
0033   ~SimpleTrackValidation() override = default;
0034
0035   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0036
0037 private:
0038   void beginJob() override;
0039   void analyze(const edm::Event&, const edm::EventSetup&) override;
0040   void endJob() override;
0041
0042   // Counters
0043   int global_rt_ = 0;
0044   int global_at_ = 0;
0045   int global_st_ = 0;
0046   int global_dt_ = 0;
0047   int global_ast_ = 0;
0048
0049   TrackingParticleSelector tpSelector;
0050   TTree* output_tree_;
0051   const std::vector<edm::InputTag> trackLabels_;
0052   std::vector<edm::EDGetTokenT<edm::View<reco::Track>>> trackTokens_;
0053   const edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator> trackAssociatorToken_;
0054   const edm::EDGetTokenT<TrackingParticleCollection> trackingParticleToken_;
0055 };
0056
0057 SimpleTrackValidation::SimpleTrackValidation(const edm::ParameterSet& iConfig)
0058     : trackLabels_(iConfig.getParameter<std::vector<edm::InputTag>>("trackLabels")),
0059       trackAssociatorToken_(consumes<reco::TrackToTrackingParticleAssociator>(
0060           iConfig.getUntrackedParameter<edm::InputTag>("trackAssociator"))),
0061       trackingParticleToken_(
0062           consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("trackingParticles"))) {
0063   for (auto& itag : trackLabels_) {
0064     trackTokens_.push_back(consumes<edm::View<reco::Track>>(itag));
0065   }
0066   tpSelector = TrackingParticleSelector(iConfig.getParameter<double>("ptMinTP"),
0067                                         iConfig.getParameter<double>("ptMaxTP"),
0068                                         iConfig.getParameter<double>("minRapidityTP"),
0069                                         iConfig.getParameter<double>("maxRapidityTP"),
0070                                         iConfig.getParameter<double>("tipTP"),
0071                                         iConfig.getParameter<double>("lipTP"),
0072                                         iConfig.getParameter<int>("minHitTP"),
0073                                         iConfig.getParameter<bool>("signalOnlyTP"),
0074                                         iConfig.getParameter<bool>("intimeOnlyTP"),
0075                                         iConfig.getParameter<bool>("chargedOnlyTP"),
0076                                         iConfig.getParameter<bool>("stableOnlyTP"),
0077                                         iConfig.getParameter<std::vector<int>>("pdgIdTP"),
0078                                         iConfig.getParameter<bool>("invertRapidityCutTP"),
0079                                         iConfig.getParameter<double>("minPhiTP"),
0080                                         iConfig.getParameter<double>("maxPhiTP"));
0081 }
0082
0083 void SimpleTrackValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0084   using namespace edm;
0085
0086   auto const& associatorByHits = iEvent.get(trackAssociatorToken_);
0087   auto TPCollectionH = iEvent.getHandle(trackingParticleToken_);
0088   TrackingParticleRefVector tpCollection;
0089
0090   for (size_t i = 0, size = TPCollectionH->size(); i < size; ++i) {
0091     auto tp = TrackingParticleRef(TPCollectionH, i);
0092     if (tpSelector(*tp)) {
0093       tpCollection.push_back(tp);
0094     }
0095   }
0096
0097   for (const auto& trackToken : trackTokens_) {
0098     edm::Handle<edm::View<reco::Track>> tracksHandle;
0099     iEvent.getByToken(trackToken, tracksHandle);
0100     const edm::View<reco::Track>& tracks = *tracksHandle;
0101
0102     edm::RefToBaseVector<reco::Track> trackRefs;
0103     for (edm::View<reco::Track>::size_type i = 0; i < tracks.size(); ++i) {
0104       trackRefs.push_back(tracks.refAt(i));
0105     }
0106
0107     reco::RecoToSimCollection recSimColl = associatorByHits.associateRecoToSim(trackRefs, tpCollection);
0108     reco::SimToRecoCollection simRecColl = associatorByHits.associateSimToReco(trackRefs, tpCollection);
0109
0110     int rt = 0;                    // number of reconstructed tracks;
0111     int at = 0;                    // number of reconstructed tracks associated to a tracking particle
0112     int ast = 0;                   // number of tracking particles associated to at least a reconstructed track.
0113     int dt = 0;                    //  number of duplicates (number of sim associated to more than one reco);
0114     int st = tpCollection.size();  // number of tracking particles;
0115
0116     for (const auto& track : trackRefs) {
0117       rt++;
0118       auto foundTP = recSimColl.find(track);
0119       if (foundTP != recSimColl.end()) {
0120         const auto& tp = foundTP->val;
0121         if (!tp.empty()) {
0122           at++;
0123         }
0124         if (simRecColl.find(tp[0].first) != simRecColl.end()) {
0125           if (simRecColl[tp[0].first].size() > 1) {
0126             dt++;
0127           }
0128         }
0129       }
0130     }
0131     for (const TrackingParticleRef& tpr : tpCollection) {
0132       auto foundTrack = simRecColl.find(tpr);
0133       if (foundTrack != simRecColl.end()) {
0134         ast++;
0135       }
0136     }
0137
0138     LogPrint("TrackValidator") << "Tag " << trackLabels_[0].label() << " Total simulated " << st
0139                                << " Associated tracks " << at << " Total reconstructed " << rt;
0140     global_rt_ += rt;
0141     global_st_ += st;
0142     global_at_ += at;
0143     global_dt_ += dt;
0144     global_ast_ += ast;
0145   }
0146 }
0147
0148 void SimpleTrackValidation::beginJob() {
0149   edm::Service<TFileService> fs;
0150   output_tree_ = fs->make<TTree>("output", "Simple Track Validation TTree");
0151
0152   output_tree_->Branch("rt", &global_rt_);
0153   output_tree_->Branch("at", &global_at_);
0154   output_tree_->Branch("st", &global_st_);
0155   output_tree_->Branch("dt", &global_dt_);
0156   output_tree_->Branch("ast", &global_ast_);
0157 }
0158
0159 void SimpleTrackValidation::endJob() { output_tree_->Fill(); }
0160
0161 void SimpleTrackValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0162   edm::ParameterSetDescription desc;
0163
0164   desc.add<std::vector<edm::InputTag>>("trackLabels", {edm::InputTag("generalTracks")});
0165   desc.add<edm::InputTag>("trackingParticles", edm::InputTag("mix", "MergedTrackTruth"));
0166   desc.add<edm::InputTag>("trackAssociator", edm::InputTag("trackingParticleRecoTrackAsssociation"));
0167
0168   // TP Selector parameters
0169   desc.add<double>("ptMinTP", 0.9);
0170   desc.add<double>("ptMaxTP", 1e100);
0171   desc.add<double>("minRapidityTP", -2.4);
0172   desc.add<double>("maxRapidityTP", 2.4);
0173   desc.add<double>("tipTP", 3.5);
0174   desc.add<double>("lipTP", 30.0);
0175   desc.add<int>("minHitTP", 0);
0176   desc.add<bool>("signalOnlyTP", true);
0177   desc.add<bool>("intimeOnlyTP", false);
0178   desc.add<bool>("chargedOnlyTP", true);
0179   desc.add<bool>("stableOnlyTP", false);
0180   desc.add<std::vector<int>>("pdgIdTP", {});
0181   desc.add<bool>("invertRapidityCutTP", false);
0182   desc.add<double>("minPhiTP", -3.2);
0183   desc.add<double>("maxPhiTP", 3.2);
0184
0185   descriptions.addWithDefaultLabel(desc);
0186 }
0187
0188 DEFINE_FWK_MODULE(SimpleTrackValidation);