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 }
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
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
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
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
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
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
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
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
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 }
0270
0271 }
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 }
0338
0339 }
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);