File indexing completed on 2023-03-31 23:02:27
0001 #include "FWCore/Framework/interface/MakerMacros.h"
0002 #include "FWCore/Framework/interface/Frameworkfwd.h"
0003 #include "FWCore/Framework/interface/EDAnalyzer.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "DataFormats/Common/interface/Handle.h"
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "FWCore/ServiceRegistry/interface/Service.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0013 #include "DataFormats/TrackReco/interface/Track.h"
0014 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0015 #include "DataFormats/VertexReco/interface/Vertex.h"
0016
0017 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0018 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0019 #include "SimDataFormats/Track/interface/SimTrack.h"
0020 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0021
0022
0023
0024 #include "RecoTracker/PixelVertexFinding/interface/PVPositionBuilder.h"
0025 #include "RecoTracker/PixelVertexFinding/interface/PVClusterComparer.h"
0026
0027 #include "MagneticField/Engine/interface/MagneticField.h"
0028 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0029
0030 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0031 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0032 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0033 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0034 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0035
0036 #include <iostream>
0037 #include <vector>
0038 #include <cmath>
0039 #include "TTree.h"
0040 #include "TFile.h"
0041 #include "TDirectory.h"
0042
0043 using namespace std;
0044
0045 class PixelVertexTest : public edm::EDAnalyzer {
0046 public:
0047 explicit PixelVertexTest(const edm::ParameterSet& conf);
0048 ~PixelVertexTest();
0049 virtual void beginJob();
0050 virtual void analyze(const edm::Event& ev, const edm::EventSetup& es);
0051 virtual void endJob();
0052
0053 private:
0054 edm::ParameterSet conf_;
0055
0056 int verbose_;
0057
0058 TTree* t_;
0059 TFile* f_;
0060 int ntrk_;
0061 static const int maxtrk_ = 1000;
0062 double pt_[maxtrk_];
0063 double z0_[maxtrk_];
0064 double errz0_[maxtrk_];
0065
0066 double theta_[maxtrk_];
0067 int nvtx_;
0068 static const int maxvtx_ = 15;
0069 double vz_[maxvtx_];
0070 double vzwt_[maxvtx_];
0071 double errvz_[maxvtx_];
0072 double errvzwt_[maxvtx_];
0073 int nvtx2_;
0074 double vz2_[maxvtx_];
0075 double trk2avg_[maxvtx_];
0076 double errvz2_[maxvtx_];
0077 int ntrk2_[maxvtx_];
0078 double sumpt2_[maxvtx_];
0079 double simx_;
0080 double simy_;
0081 double simz_;
0082 double vxkal_[maxvtx_];
0083 double vykal_[maxvtx_];
0084 double vzkal_[maxvtx_];
0085 };
0086
0087 PixelVertexTest::PixelVertexTest(const edm::ParameterSet& conf) : conf_(conf), t_(0), f_(0) {
0088 edm::LogInfo("PixelVertexTest") << " CTOR";
0089 }
0090
0091 PixelVertexTest::~PixelVertexTest() {
0092 edm::LogInfo("PixelVertexTest") << " DTOR";
0093 delete f_;
0094
0095 }
0096
0097 void PixelVertexTest::beginJob() {
0098
0099 verbose_ = conf_.getUntrackedParameter<unsigned int>("Verbosity", 0);
0100
0101
0102 std::string file = conf_.getUntrackedParameter<std::string>("OutputTree", "mytree.root");
0103 const char* cwd = gDirectory->GetPath();
0104 f_ = new TFile(file.c_str(), "RECREATE");
0105 t_ = new TTree("t", "Pixel Vertex Testing");
0106 t_->Branch("nvtx", &nvtx_, "nvtx/I");
0107 t_->Branch("vz", vz_, "vz[nvtx]/D");
0108 t_->Branch("errvz", errvz_, "errvz[nvtx]/D");
0109 t_->Branch("vzwt", vzwt_, "vzwt[nvtx]/D");
0110 t_->Branch("errvzwt", errvzwt_, "errvzwt[nvtx]/D");
0111 t_->Branch("nvtx2", &nvtx2_, "nvtx2/I");
0112 t_->Branch("vz2", vz2_, "vz2[nvtx2]/D");
0113 t_->Branch("trk2avg", trk2avg_, "trk2avg[nvtx2]/D");
0114 t_->Branch("errvz2", errvz2_, "errvz2[nvtx2]/D");
0115 t_->Branch("ntrk2", ntrk2_, "ntrk2[nvtx2]/I");
0116 t_->Branch("sumpt2", sumpt2_, "sumpt2[nvtx2]/D");
0117 t_->Branch("ntrk", &ntrk_, "ntrk/I");
0118 t_->Branch("pt", pt_, "pt[ntrk]/D");
0119 t_->Branch("z0", z0_, "z0[ntrk]/D");
0120 t_->Branch("errz0", errz0_, "errz0[ntrk]/D");
0121
0122 t_->Branch("theta", theta_, "theta[ntrk]/D");
0123 t_->Branch("simx", &simx_, "simx/D");
0124 t_->Branch("simy", &simy_, "simy/D");
0125 t_->Branch("simz", &simz_, "simz/D");
0126 t_->Branch("vxkal", vxkal_, "vxkal[nvtx2]/D");
0127 t_->Branch("vykal", vykal_, "vykal[nvtx2]/D");
0128 t_->Branch("vzkal", vzkal_, "vzkal[nvtx2]/D");
0129 gDirectory->cd(cwd);
0130 }
0131
0132 void PixelVertexTest::analyze(const edm::Event& ev, const edm::EventSetup& es) {
0133 cout << "*** PixelVertexTest, analyze event: " << ev.id() << endl;
0134
0135 edm::ESHandle<MagneticField> field;
0136 es.get<IdealMagneticFieldRecord>().get(field);
0137
0138 edm::InputTag simG4 = conf_.getParameter<edm::InputTag>("simG4");
0139 edm::Handle<edm::SimVertexContainer> simVtcs;
0140 ev.getByLabel(simG4, simVtcs);
0141 if (verbose_ > 0) {
0142 cout << "simulated vertices: " << simVtcs->size() << std::endl;
0143 }
0144
0145
0146
0147
0148 simx_ = (simVtcs->size() > 0) ? (*simVtcs)[0].position().x() : -9999.0;
0149 simy_ = (simVtcs->size() > 0) ? (*simVtcs)[0].position().y() : -9999.0;
0150 simz_ = (simVtcs->size() > 0) ? (*simVtcs)[0].position().z() : -9999.0;
0151 if (verbose_ > 1) {
0152 for (int i = 0; i < simVtcs->size(); i++) {
0153 std::cout << (*simVtcs)[i].parentIndex() << ": " << (*simVtcs)[i].position().x() << ", "
0154 << (*simVtcs)[i].position().y() << ", " << (*simVtcs)[i].position().z() << "; ";
0155 }
0156 std::cout << "\n" << std::endl;
0157 }
0158
0159 edm::Handle<reco::TrackCollection> trackCollection;
0160 std::string trackCollName = conf_.getParameter<std::string>("TrackCollection");
0161 ev.getByLabel(trackCollName, trackCollection);
0162 const reco::TrackCollection tracks = *(trackCollection.product());
0163
0164 reco::TrackRefVector trks;
0165
0166 if (verbose_ > 0) {
0167 std::cout << *(trackCollection.provenance()) << std::endl;
0168 cout << "Reconstructed " << tracks.size() << " tracks" << std::endl;
0169 }
0170 ntrk_ = 0;
0171 for (unsigned int i = 0; i < tracks.size(); i++) {
0172 if (verbose_ > 0) {
0173 cout << "\tmomentum: " << tracks[i].momentum() << "\tPT: " << tracks[i].pt() << endl;
0174 cout << "\tvertex: " << tracks[i].vertex() << "\tZ0: " << tracks[i].dz() << " +- " << tracks[i].dzError() << endl;
0175 cout << "\tcharge: " << tracks[i].charge() << endl;
0176 }
0177 trks.push_back(reco::TrackRef(trackCollection, i));
0178
0179 if (ntrk_ < maxtrk_) {
0180 pt_[ntrk_] = tracks[i].pt();
0181 z0_[ntrk_] = tracks[i].dz();
0182
0183 errz0_[ntrk_] = tracks[i].dzError();
0184
0185 theta_[ntrk_] = tracks[i].theta();
0186 ntrk_++;
0187 }
0188 if (verbose_ > 0)
0189 cout << "------------------------------------------------" << endl;
0190 }
0191 PVPositionBuilder pos;
0192 nvtx_ = 0;
0193 vz_[nvtx_] = pos.average(trks).value();
0194 errvz_[nvtx_] = pos.average(trks).error();
0195 vzwt_[nvtx_] = pos.wtAverage(trks).value();
0196 errvzwt_[nvtx_] = pos.wtAverage(trks).error();
0197 nvtx_++;
0198 if (verbose_ > 0) {
0199 std::cout << "The average z-position of these tracks is " << vz_[0] << " +- " << errvz_[0] << std::endl;
0200 std::cout << "The weighted average z-position of these tracks is " << vzwt_[0] << " +- " << errvzwt_[0]
0201 << std::endl;
0202 }
0203
0204
0205 edm::Handle<reco::VertexCollection> vertexCollection;
0206 ev.getByLabel("pixelVertices", vertexCollection);
0207 const reco::VertexCollection vertexes = *(vertexCollection.product());
0208 if (verbose_ > 0) {
0209 std::cout << *(vertexCollection.provenance()) << std::endl;
0210 cout << "Reconstructed " << vertexes.size() << " vertexes" << std::endl;
0211 }
0212 nvtx2_ = vertexes.size();
0213 PVClusterComparer vcompare;
0214 for (int i = 0; i < nvtx2_ && i < maxvtx_; i++) {
0215 vz2_[i] = vertexes[i].z();
0216 errvz2_[i] = std::sqrt(vertexes[i].error(2, 2));
0217 ntrk2_[i] = vertexes[i].tracksSize();
0218 sumpt2_[i] = vcompare.pTSquaredSum(vertexes[i]);
0219
0220
0221 while (!trks.empty())
0222 trks.erase(trks.begin());
0223 for (reco::track_iterator j = vertexes[i].tracks_begin(); j != vertexes[i].tracks_end(); ++j)
0224 trks.push_back(*j);
0225 trk2avg_[i] = pos.wtAverage(trks).value();
0226 }
0227
0228
0229 if (nvtx2_ > 0) {
0230 vector<reco::TransientTrack> t_tks;
0231 for (int i = 0; i < nvtx2_ && i < maxvtx_; i++) {
0232 t_tks.clear();
0233 for (reco::track_iterator j = vertexes[i].tracks_begin(); j != vertexes[i].tracks_end(); ++j) {
0234 t_tks.push_back(reco::TransientTrack(**j, field.product()));
0235 }
0236 KalmanVertexFitter kvf;
0237
0238 TransientVertex tv = kvf.vertex(t_tks);
0239 if (verbose_ > 0)
0240 std::cout << "Kalman Position: " << reco::Vertex::Point(tv.position()) << std::endl;
0241 vxkal_[i] = tv.position().x();
0242 vykal_[i] = tv.position().y();
0243 vzkal_[i] = tv.position().z();
0244 }
0245 }
0246
0247
0248 t_->Fill();
0249 }
0250
0251 void PixelVertexTest::endJob() {
0252 if (t_)
0253 t_->Print();
0254 if (f_) {
0255 f_->Print();
0256 f_->Write();
0257 }
0258 }
0259
0260 DEFINE_FWK_MODULE(PixelVertexTest);