Back to home page

Project CMSSW displayed by LXR
 
 

    

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 //#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
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   // How noisy should I be
0056   int verbose_;
0057   // Tree of simple vars for testing resolution eff etc
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   //  double tanl_[maxtrk_];
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   //  delete t_;
0095 }
0096 
0097 void PixelVertexTest::beginJob() {
0098   // How noisy?
0099   verbose_ = conf_.getUntrackedParameter<unsigned int>("Verbosity", 0);
0100 
0101   // Make my little tree
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   //  t_->Branch("tanl",tanl_,"tanl[ntrk]/D");
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   //  simx_ = (simVtcs->size() > 0) ? (*simVtcs)[0].position().x()/10 : -9999.0;
0145   //  simy_ = (simVtcs->size() > 0) ? (*simVtcs)[0].position().y()/10 : -9999.0;
0146   //  simz_ = (simVtcs->size() > 0) ? (*simVtcs)[0].position().z()/10 : -9999.0;
0147   // No longer need to convert from mm as of version 1_2_0
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     // Fill ntuple vars
0179     if (ntrk_ < maxtrk_) {
0180       pt_[ntrk_] = tracks[i].pt();
0181       z0_[ntrk_] = tracks[i].dz();
0182       //      errz0_[ntrk_] = std::sqrt( tracks[i].covariance(3,3) );
0183       errz0_[ntrk_] = tracks[i].dzError();
0184       //      tanl_[ntrk_] = tracks[i].tanDip();
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   // NOW let's see if my vertex producer did a darn thing...
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     // Now calculate my own average position by hand to cross check conversion process
0220     //    trks.clear(); // not yet implemented
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   // Now let's send off our tracks to the Kalman fitter to see what great things happen...
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       //      TransientVertex tv = kvf.vertex(t_tks,GlobalPoint(vertexes[i].x(),vertexes[i].y(),vertexes[i].z()));
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   // Finally, fill the tree with the above values
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!