Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:31:08

0001 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 
0006 #include "MagneticField/Engine/interface/MagneticField.h"
0007 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0008 
0009 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0010 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0011 #include "FWCore/ServiceRegistry/interface/Service.h"
0012 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0013 #include "SimDataFormats/Associations/interface/VertexToTrackingVertexAssociator.h"
0014 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0015 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0016 
0017 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0018 #include "DataFormats/Math/interface/Vector.h"
0019 #include "DataFormats/Math/interface/LorentzVector.h"
0020 #include <Math/GenVector/PxPyPzE4D.h>
0021 #include <Math/GenVector/PxPyPzM4D.h>
0022 
0023 #include "TFile.h"
0024 #include "TH1F.h"
0025 #include "TMath.h"
0026 #include "TROOT.h"
0027 #include "TTree.h"
0028 
0029 #include <cmath>
0030 #include <iostream>
0031 #include <map>
0032 #include <memory>
0033 #include <set>
0034 #include <string>
0035 #include <vector>
0036 
0037 namespace reco {
0038   class TrackToTrackingParticleAssociator;
0039   class VertexToTrackingVertexAssociator;
0040 }  // namespace reco
0041 
0042 class testVertexAssociator : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0043 public:
0044   testVertexAssociator(const edm::ParameterSet &conf);
0045   ~testVertexAssociator() override = default;
0046   void beginJob() override;
0047   void endJob() override;
0048   void analyze(const edm::Event &, const edm::EventSetup &) override;
0049 
0050 private:
0051   const reco::TrackToTrackingParticleAssociator *associatorByChi2;
0052   const reco::TrackToTrackingParticleAssociator *associatorByHits;
0053   const reco::VertexToTrackingVertexAssociator *associatorByTracks;
0054 
0055   const edm::EDGetTokenT<reco::VertexToTrackingVertexAssociator> associatorByTracksToken;
0056   const edm::InputTag vertexCollection_;
0057   const edm::EDGetTokenT<TrackingVertexCollection> tokenTV_;
0058   const edm::EDGetTokenT<edm::View<reco::Vertex>> tokenVtx_;
0059   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tokenMF_;
0060 
0061   int n_event_;
0062   int n_rs_vertices_;
0063   int n_rs_vtxassocs_;
0064   int n_sr_vertices_;
0065   int n_sr_vtxassocs_;
0066 
0067   //--------- RecoToSim Histos -----
0068 
0069   TH1F *rs_resx;
0070   TH1F *rs_resy;
0071   TH1F *rs_resz;
0072   TH1F *rs_pullx;
0073   TH1F *rs_pully;
0074   TH1F *rs_pullz;
0075   TH1F *rs_dist;
0076   TH1F *rs_simz;
0077   TH1F *rs_recz;
0078   TH1F *rs_nrectrk;
0079   TH1F *rs_nsimtrk;
0080   TH1F *rs_qual;
0081   TH1F *rs_chi2norm;
0082   TH1F *rs_chi2prob;
0083 
0084   //--------- SimToReco Histos -----
0085 
0086   TH1F *sr_resx;
0087   TH1F *sr_resy;
0088   TH1F *sr_resz;
0089   TH1F *sr_pullx;
0090   TH1F *sr_pully;
0091   TH1F *sr_pullz;
0092   TH1F *sr_dist;
0093   TH1F *sr_simz;
0094   TH1F *sr_recz;
0095   TH1F *sr_nrectrk;
0096   TH1F *sr_nsimtrk;
0097   TH1F *sr_qual;
0098   TH1F *sr_chi2norm;
0099   TH1F *sr_chi2prob;
0100 };
0101 
0102 // class TrackAssociator;
0103 class TrackAssociatorByHits;
0104 class TrackerHitAssociator;
0105 
0106 testVertexAssociator::testVertexAssociator(edm::ParameterSet const &conf)
0107     : associatorByTracksToken(consumes<reco::VertexToTrackingVertexAssociator>(
0108           conf.getUntrackedParameter<edm::InputTag>("vertexAssociation"))),
0109       vertexCollection_(conf.getUntrackedParameter<edm::InputTag>("vertexCollection")),
0110       tokenTV_(consumes<TrackingVertexCollection>(edm::InputTag("mix", "MergedTrackTruth"))),
0111       tokenVtx_(consumes<edm::View<reco::Vertex>>(vertexCollection_)),
0112       tokenMF_(esConsumes<MagneticField, IdealMagneticFieldRecord>()) {
0113   usesResource(TFileService::kSharedResource);
0114 
0115   n_event_ = 0;
0116   n_rs_vertices_ = 0;
0117   n_rs_vtxassocs_ = 0;
0118   n_sr_vertices_ = 0;
0119   n_sr_vtxassocs_ = 0;
0120 }
0121 
0122 void testVertexAssociator::beginJob() {
0123   edm::Service<TFileService> fs;
0124 
0125   // RecoToSim Histos
0126 
0127   rs_dist = fs->make<TH1F>("rs_dist", "r Miss Distance (cm)", 100, 0, 0.1);
0128   rs_recz = fs->make<TH1F>("rs_recz", "z, Reconstructed Vertex (cm)", 200, -25.0, 25.0);
0129   rs_simz = fs->make<TH1F>("rs_simz", "z, Simulated Vertex (cm)", 200, -25.0, 25.0);
0130   rs_nsimtrk = fs->make<TH1F>("rs_nsimtrk", "# of tracks, Simulated", 501, -0.5, 500.5);
0131   rs_nrectrk = fs->make<TH1F>("rs_nrectrk", "# of tracks, Reconstructed", 501, -0.5, 500.5);
0132   rs_qual = fs->make<TH1F>("rs_qual", "Quality of Match", 51, -0.01, 1.01);
0133   rs_chi2norm = fs->make<TH1F>("rs_chi2norm", "Vertex Normalized Chi2", 100, 0, 10.);
0134   rs_chi2prob = fs->make<TH1F>("rs_chi2prob", "Vertex Chi2 Probability", 100, 0, 1.);
0135   rs_resx = fs->make<TH1F>("rs_resx", "rec - sim Distance (cm)", 100, -0.05, 0.05);
0136   rs_resy = fs->make<TH1F>("rs_resy", "rec - sim Distance (cm)", 100, -0.05, 0.05);
0137   rs_resz = fs->make<TH1F>("rs_resz", "rec - sim Distance (cm)", 100, -0.05, 0.05);
0138   rs_pullx = fs->make<TH1F>("rs_pullx", "(rec - sim)/err_rec ", 100, -10., 10.);
0139   rs_pully = fs->make<TH1F>("rs_pully", "(rec - sim)/err_rec ", 100, -10., 10.);
0140   rs_pullz = fs->make<TH1F>("rs_pullz", "(rec - sim)/err_rec ", 100, -10., 10.);
0141 
0142   // SimToReco Histos
0143 
0144   sr_dist = fs->make<TH1F>("sr_dist", "r Miss Distance (cm)", 100, 0, 0.1);
0145   sr_recz = fs->make<TH1F>("sr_recz", "z, Reconstructed Vertex (cm)", 200, -25.0, 25.0);
0146   sr_simz = fs->make<TH1F>("sr_simz", "z, Simulated Vertex (cm)", 200, -25.0, 25.0);
0147   sr_nsimtrk = fs->make<TH1F>("sr_nsimtrk", "# of tracks, Simulated", 501, -0.5, 500.5);
0148   sr_nrectrk = fs->make<TH1F>("sr_nrectrk", "# of tracks, Reconstructed", 501, -0.5, 500.5);
0149   sr_qual = fs->make<TH1F>("sr_qual", "Quality of Match", 51, -0.01, 1.01);
0150   sr_chi2norm = fs->make<TH1F>("sr_chi2norm", "Vertex Normalized Chi2", 100, 0, 10.);
0151   sr_chi2prob = fs->make<TH1F>("sr_chi2prob", "Vertex Chi2 Probability", 100, 0, 1.);
0152   sr_resx = fs->make<TH1F>("sr_resx", "rec - sim Distance (cm)", 100, -0.05, 0.05);
0153   sr_resy = fs->make<TH1F>("sr_resy", "rec - sim Distance (cm)", 100, -0.05, 0.05);
0154   sr_resz = fs->make<TH1F>("sr_resz", "rec - sim Distance (cm)", 100, -0.05, 0.05);
0155   sr_pullx = fs->make<TH1F>("sr_pullx", "(rec - sim)/err_rec ", 100, -10., 10.);
0156   sr_pully = fs->make<TH1F>("sr_pully", "(rec - sim)/err_rec ", 100, -10., 10.);
0157   sr_pullz = fs->make<TH1F>("sr_pullz", "(rec - sim)/err_rec ", 100, -10., 10.);
0158 }
0159 
0160 void testVertexAssociator::endJob() {
0161   std::cout << std::endl;
0162   std::cout << " ====== Total Number of analyzed events: " << n_event_ << " ====== " << std::endl;
0163   std::cout << " ====== Total Number of R2S vertices:    " << n_rs_vertices_ << " ====== " << std::endl;
0164   std::cout << " ====== Total Number of R2S vtx assocs:  " << n_rs_vtxassocs_ << " ====== " << std::endl;
0165   std::cout << " ====== Total Number of S2R vertices:    " << n_sr_vertices_ << " ====== " << std::endl;
0166   std::cout << " ====== Total Number of S2R vtx assocs:  " << n_sr_vtxassocs_ << " ====== " << std::endl;
0167 }
0168 
0169 void testVertexAssociator::analyze(const edm::Event &event, const edm::EventSetup &setup) {
0170   //const auto &theMF = setup.getHandle(tokenMF_);
0171 
0172   const edm::Handle<reco::VertexToTrackingVertexAssociator> &theTracksAssociator =
0173       event.getHandle(associatorByTracksToken);
0174   associatorByTracks = theTracksAssociator.product();
0175 
0176   ++n_event_;
0177 
0178   std::cout << "*** Analyzing " << event.id() << " n_event = " << n_event_ << std::endl << std::endl;
0179 
0180   const auto &TVCollection = event.getHandle(tokenTV_);
0181   const TrackingVertexCollection tVC = *(TVCollection.product());
0182 
0183   // Vertex Collection
0184   const edm::Handle<edm::View<reco::Vertex>> &vertexCollection = event.getHandle(tokenVtx_);
0185   const edm::View<reco::Vertex> vC = *(vertexCollection.product());
0186 
0187   std::cout << std::endl;
0188   std::cout << "                      ****************** Before Assocs "
0189                "****************** "
0190             << std::endl
0191             << std::endl;
0192 
0193   std::cout << "vertexCollection.size() = " << vC.size() << std::endl;
0194   std::cout << "TVCollection.size()     = " << tVC.size() << std::endl;
0195 
0196   std::cout << std::endl;
0197   std::cout << "                      ****************** Reco To Sim "
0198                "****************** "
0199             << std::endl
0200             << std::endl;
0201 
0202   // std::cout << "-- Associator by hits --" << std::endl;
0203   reco::VertexRecoToSimCollection r2sVertices = associatorByTracks->associateRecoToSim(vertexCollection, TVCollection);
0204 
0205   reco::VertexSimToRecoCollection s2rVertices = associatorByTracks->associateSimToReco(vertexCollection, TVCollection);
0206 
0207   std::cout << std::endl;
0208   std::cout << "VertexRecoToSim size           = " << r2sVertices.size()
0209             << " ; VertexSimToReco size           = " << r2sVertices.size() << " " << std::endl;
0210 
0211   std::cout << std::endl << " [testVertexAssociator] Analyzing Reco To Sim" << std::endl;
0212 
0213   int cont_recvR2S = 0;
0214 
0215   for (reco::VertexRecoToSimCollection::const_iterator iR2S = r2sVertices.begin(); iR2S != r2sVertices.end();
0216        ++iR2S, ++cont_recvR2S) {
0217     ++n_rs_vertices_;
0218 
0219     reco::VertexBaseRef recVertex = iR2S->key;
0220     math::XYZPoint recPos = recVertex->position();
0221 
0222     double nrectrk = recVertex->tracksSize();
0223 
0224     std::vector<std::pair<TrackingVertexRef, double>> simVertices = iR2S->val;
0225 
0226     int cont_simvR2S = 0;
0227     for (std::vector<std::pair<TrackingVertexRef, double>>::const_iterator iMatch = simVertices.begin();
0228          iMatch != simVertices.end();
0229          ++iMatch, ++cont_simvR2S) {
0230       TrackingVertexRef simVertex = iMatch->first;
0231       math::XYZTLorentzVectorD simVec = (iMatch->first)->position();
0232       math::XYZPoint simPos = math::XYZPoint(simVec.x(), simVec.y(), simVec.z());
0233 
0234       ++n_rs_vtxassocs_;
0235 
0236       std::cout << "rec vertex " << cont_recvR2S << " has associated sim vertex " << cont_simvR2S << std::endl;
0237 
0238       double nsimtrk = simVertex->daughterTracks().size();
0239       double qual = iMatch->second;
0240 
0241       double chi2norm = recVertex->normalizedChi2();
0242       double chi2prob = ChiSquaredProbability(recVertex->chi2(), recVertex->ndof());
0243 
0244       double resx = recVertex->x() - simVertex->position().x();
0245       double resy = recVertex->y() - simVertex->position().y();
0246       double resz = recVertex->z() - simVertex->position().z();
0247       double pullx = (recVertex->x() - simVertex->position().x()) / recVertex->xError();
0248       double pully = (recVertex->y() - simVertex->position().y()) / recVertex->yError();
0249       double pullz = (recVertex->z() - simVertex->position().z()) / recVertex->zError();
0250       double dist = sqrt(resx * resx + resy * resy + resz * resz);
0251 
0252       std::cout << "            R2S: recPos = " << recPos << " ; simPos = " << simPos << std::endl;
0253 
0254       rs_resx->Fill(resx);
0255       rs_resy->Fill(resy);
0256       rs_resz->Fill(resz);
0257       rs_pullx->Fill(pullx);
0258       rs_pully->Fill(pully);
0259       rs_pullz->Fill(pullz);
0260       rs_dist->Fill(dist);
0261       rs_simz->Fill(simPos.Z());
0262       rs_recz->Fill(recPos.Z());
0263       rs_nsimtrk->Fill(nsimtrk);
0264       rs_nrectrk->Fill(nrectrk);
0265       rs_qual->Fill(qual);
0266       rs_chi2norm->Fill(chi2norm);
0267       rs_chi2prob->Fill(chi2prob);
0268 
0269     }  // end simVertices
0270 
0271   }  // end iR2S
0272 
0273   std::cout << std::endl
0274             << "                      ****************** Sim To Reco "
0275                "****************** "
0276             << std::endl
0277             << std::endl;
0278 
0279   std::cout << std::endl << "[testVertexAssociator] Analyzing Sim To Reco" << std::endl;
0280 
0281   int cont_simvS2R = 0;
0282   for (reco::VertexSimToRecoCollection::const_iterator iS2R = s2rVertices.begin(); iS2R != s2rVertices.end();
0283        ++iS2R, ++cont_simvS2R) {
0284     ++n_sr_vertices_;
0285 
0286     TrackingVertexRef simVertex = iS2R->key;
0287     math::XYZTLorentzVectorD simVec = simVertex->position();
0288     math::XYZPoint simPos = math::XYZPoint(simVec.x(), simVec.y(), simVec.z());
0289 
0290     double nsimtrk = simVertex->daughterTracks().size();
0291 
0292     std::vector<std::pair<reco::VertexBaseRef, double>> recoVertices = iS2R->val;
0293 
0294     int cont_recvS2R = 0;
0295 
0296     for (std::vector<std::pair<reco::VertexBaseRef, double>>::const_iterator iMatch = recoVertices.begin();
0297          iMatch != recoVertices.end();
0298          ++iMatch, ++cont_recvS2R) {
0299       reco::VertexBaseRef recVertex = iMatch->first;
0300       math::XYZPoint recPos = recVertex->position();
0301 
0302       ++n_sr_vtxassocs_;
0303 
0304       std::cout << "sim vertex " << cont_simvS2R << " has associated rec vertex " << cont_recvS2R << std::endl;
0305 
0306       double nrectrk = recVertex->tracksSize();
0307       double qual = iMatch->second;
0308 
0309       double chi2norm = recVertex->normalizedChi2();
0310       double chi2prob = ChiSquaredProbability(recVertex->chi2(), recVertex->ndof());
0311 
0312       double resx = recVertex->x() - simVertex->position().x();
0313       double resy = recVertex->y() - simVertex->position().y();
0314       double resz = recVertex->z() - simVertex->position().z();
0315       double pullx = (recVertex->x() - simVertex->position().x()) / recVertex->xError();
0316       double pully = (recVertex->y() - simVertex->position().y()) / recVertex->yError();
0317       double pullz = (recVertex->z() - simVertex->position().z()) / recVertex->zError();
0318       double dist = sqrt(resx * resx + resy * resy + resz * resz);
0319 
0320       std::cout << "            S2R: simPos = " << simPos << " ; recPos = " << recPos << std::endl;
0321 
0322       sr_resx->Fill(resx);
0323       sr_resy->Fill(resy);
0324       sr_resz->Fill(resz);
0325       sr_pullx->Fill(pullx);
0326       sr_pully->Fill(pully);
0327       sr_pullz->Fill(pullz);
0328       sr_dist->Fill(dist);
0329       sr_simz->Fill(simPos.Z());
0330       sr_recz->Fill(recPos.Z());
0331       sr_nsimtrk->Fill(nsimtrk);
0332       sr_nrectrk->Fill(nrectrk);
0333       sr_qual->Fill(qual);
0334       sr_chi2norm->Fill(chi2norm);
0335       sr_chi2prob->Fill(chi2prob);
0336 
0337     }  // end recoVertices
0338 
0339   }  // end iS2R
0340 
0341   std::cout << std::endl;
0342 }
0343 
0344 #include "FWCore/Framework/interface/MakerMacros.h"
0345 #include "FWCore/PluginManager/interface/ModuleDef.h"
0346 DEFINE_FWK_MODULE(testVertexAssociator);