Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:26:52

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 "RecoPixelVertexing/PixelVertexFinding/interface/PVPositionBuilder.h"
0025 #include "RecoPixelVertexing/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 private:
0053   edm::ParameterSet conf_; 
0054   // How noisy should I be
0055   int verbose_;
0056   // Tree of simple vars for testing resolution eff etc
0057   TTree *t_;
0058   TFile *f_;
0059   int ntrk_;
0060   static const int maxtrk_=1000;
0061   double pt_[maxtrk_];
0062   double z0_[maxtrk_];
0063   double errz0_[maxtrk_];
0064   //  double tanl_[maxtrk_];
0065   double theta_[maxtrk_];
0066   int nvtx_;
0067   static const int maxvtx_=15;
0068   double vz_[maxvtx_];
0069   double vzwt_[maxvtx_];
0070   double errvz_[maxvtx_];
0071   double errvzwt_[maxvtx_];
0072   int nvtx2_;
0073   double vz2_[maxvtx_];
0074   double trk2avg_[maxvtx_];
0075   double errvz2_[maxvtx_];
0076   int ntrk2_[maxvtx_];
0077   double sumpt2_[maxvtx_];
0078   double simx_;
0079   double simy_;
0080   double simz_;
0081   double vxkal_[maxvtx_];
0082   double vykal_[maxvtx_];
0083   double vzkal_[maxvtx_];
0084 };
0085 
0086 PixelVertexTest::PixelVertexTest(const edm::ParameterSet& conf)
0087   : conf_(conf),t_(0),f_(0)
0088 {
0089   edm::LogInfo("PixelVertexTest")<<" CTOR";
0090 }
0091 
0092 PixelVertexTest::~PixelVertexTest()
0093 {
0094   edm::LogInfo("PixelVertexTest")<<" DTOR";
0095   delete f_;
0096   //  delete t_;
0097 }
0098 
0099 void PixelVertexTest::beginJob() {
0100   // How noisy?
0101   verbose_ = conf_.getUntrackedParameter<unsigned int>("Verbosity",0);
0102 
0103   // Make my little tree
0104   std::string file = conf_.getUntrackedParameter<std::string>("OutputTree","mytree.root");
0105   const char* cwd= gDirectory->GetPath();
0106   f_ = new TFile(file.c_str(),"RECREATE");
0107   t_ = new TTree("t","Pixel Vertex Testing");
0108   t_->Branch("nvtx",&nvtx_,"nvtx/I");
0109   t_->Branch("vz",vz_,"vz[nvtx]/D");
0110   t_->Branch("errvz",errvz_,"errvz[nvtx]/D");
0111   t_->Branch("vzwt",vzwt_,"vzwt[nvtx]/D");
0112   t_->Branch("errvzwt",errvzwt_,"errvzwt[nvtx]/D");
0113   t_->Branch("nvtx2",&nvtx2_,"nvtx2/I");
0114   t_->Branch("vz2",vz2_,"vz2[nvtx2]/D");
0115   t_->Branch("trk2avg",trk2avg_,"trk2avg[nvtx2]/D");
0116   t_->Branch("errvz2",errvz2_,"errvz2[nvtx2]/D");
0117   t_->Branch("ntrk2",ntrk2_,"ntrk2[nvtx2]/I");
0118   t_->Branch("sumpt2",sumpt2_,"sumpt2[nvtx2]/D");
0119   t_->Branch("ntrk",&ntrk_,"ntrk/I");
0120   t_->Branch("pt",pt_,"pt[ntrk]/D");
0121   t_->Branch("z0",z0_,"z0[ntrk]/D");
0122   t_->Branch("errz0",errz0_,"errz0[ntrk]/D");
0123   //  t_->Branch("tanl",tanl_,"tanl[ntrk]/D");
0124   t_->Branch("theta",theta_,"theta[ntrk]/D");
0125   t_->Branch("simx",&simx_,"simx/D");
0126   t_->Branch("simy",&simy_,"simy/D");
0127   t_->Branch("simz",&simz_,"simz/D");
0128   t_->Branch("vxkal",vxkal_,"vxkal[nvtx2]/D");
0129   t_->Branch("vykal",vykal_,"vykal[nvtx2]/D");
0130   t_->Branch("vzkal",vzkal_,"vzkal[nvtx2]/D");
0131   gDirectory->cd(cwd);
0132 }
0133 
0134 void PixelVertexTest::analyze(
0135     const edm::Event& ev, const edm::EventSetup& es)
0136 {
0137   cout <<"*** PixelVertexTest, analyze event: " << ev.id() << endl;
0138 
0139   edm::ESHandle<MagneticField> field;
0140   es.get<IdealMagneticFieldRecord>().get(field);
0141 
0142   edm::InputTag simG4 = conf_.getParameter<edm::InputTag>( "simG4" );
0143   edm::Handle<edm::SimVertexContainer> simVtcs;
0144   ev.getByLabel( simG4, simVtcs);
0145   if (verbose_ > 0) {
0146     cout << "simulated vertices: "<< simVtcs->size() << std::endl;
0147   }
0148 //  simx_ = (simVtcs->size() > 0) ? (*simVtcs)[0].position().x()/10 : -9999.0;
0149 //  simy_ = (simVtcs->size() > 0) ? (*simVtcs)[0].position().y()/10 : -9999.0;
0150 //  simz_ = (simVtcs->size() > 0) ? (*simVtcs)[0].position().z()/10 : -9999.0;
0151 // No longer need to convert from mm as of version 1_2_0
0152   simx_ = (simVtcs->size() > 0) ? (*simVtcs)[0].position().x() : -9999.0;
0153   simy_ = (simVtcs->size() > 0) ? (*simVtcs)[0].position().y() : -9999.0;
0154   simz_ = (simVtcs->size() > 0) ? (*simVtcs)[0].position().z() : -9999.0;
0155   if (verbose_ > 1) {
0156     for (int i=0; i<simVtcs->size(); i++) {
0157       std::cout << (*simVtcs)[i].parentIndex() << ": " << (*simVtcs)[i].position().x() << ", " << (*simVtcs)[i].position().y() << ", " << (*simVtcs)[i].position().z() << ";  ";
0158     }
0159     std::cout << "\n" << std::endl;
0160   }
0161 
0162 
0163   edm::Handle<reco::TrackCollection> trackCollection;
0164   std::string trackCollName = conf_.getParameter<std::string>("TrackCollection");
0165   ev.getByLabel(trackCollName,trackCollection);
0166   const reco::TrackCollection tracks = *(trackCollection.product());
0167 
0168   reco::TrackRefVector trks;
0169 
0170   if (verbose_ > 0) {
0171     std::cout << *(trackCollection.provenance()) << std::endl;
0172     cout << "Reconstructed "<< tracks.size() << " tracks" << std::endl;
0173   }
0174   ntrk_=0;
0175   for (unsigned int i=0; i<tracks.size(); i++) {
0176     if (verbose_ > 0) {
0177       cout << "\tmomentum: " << tracks[i].momentum()
0178        << "\tPT: " << tracks[i].pt()<< endl;
0179       cout << "\tvertex: " << tracks[i].vertex()
0180        << "\tZ0: " << tracks[i].dz() << " +- " << tracks[i].dzError() << endl;
0181       cout << "\tcharge: " << tracks[i].charge()<< endl;
0182     }
0183     trks.push_back( reco::TrackRef(trackCollection, i) );
0184     // Fill ntuple vars
0185     if (ntrk_ < maxtrk_) {
0186       pt_[ntrk_] = tracks[i].pt();
0187       z0_[ntrk_] = tracks[i].dz();
0188       //      errz0_[ntrk_] = std::sqrt( tracks[i].covariance(3,3) );
0189       errz0_[ntrk_] = tracks[i].dzError();
0190       //      tanl_[ntrk_] = tracks[i].tanDip();
0191       theta_[ntrk_] = tracks[i].theta();
0192       ntrk_++;
0193     }
0194     if (verbose_ > 0) cout <<"------------------------------------------------"<<endl;
0195   }
0196   PVPositionBuilder pos;  
0197   nvtx_ = 0;
0198   vz_[nvtx_] = pos.average(trks).value();
0199   errvz_[nvtx_] = pos.average(trks).error();
0200   vzwt_[nvtx_] = pos.wtAverage(trks).value();
0201   errvzwt_[nvtx_] = pos.wtAverage(trks).error();
0202   nvtx_++;
0203   if (verbose_ > 0) {
0204     std::cout << "The average z-position of these tracks is " << vz_[0] << " +- " << errvz_[0] << std::endl;
0205     std::cout << "The weighted average z-position of these tracks is " << vzwt_[0] << " +- " << errvzwt_[0] << std::endl;
0206   }
0207 
0208   // NOW let's see if my vertex producer did a darn thing...
0209   edm::Handle<reco::VertexCollection> vertexCollection;
0210   ev.getByLabel("pixelVertices",vertexCollection);
0211   const reco::VertexCollection vertexes = *(vertexCollection.product());
0212   if (verbose_ > 0) {
0213     std::cout << *(vertexCollection.provenance()) << std::endl;
0214     cout << "Reconstructed "<< vertexes.size() << " vertexes" << std::endl;
0215   }
0216   nvtx2_ = vertexes.size();
0217   PVClusterComparer vcompare;
0218   for (int i=0; i<nvtx2_ && i<maxvtx_; i++) {
0219     vz2_[i] = vertexes[i].z();
0220     errvz2_[i] = std::sqrt(vertexes[i].error(2,2));
0221     ntrk2_[i] = vertexes[i].tracksSize();
0222     sumpt2_[i] = vcompare.pTSquaredSum(vertexes[i]);
0223     // Now calculate my own average position by hand to cross check conversion process
0224     //    trks.clear(); // not yet implemented
0225     while (! trks.empty()) trks.erase( trks.begin() );
0226     for (reco::track_iterator j=vertexes[i].tracks_begin(); j!=vertexes[i].tracks_end(); ++j) 
0227       trks.push_back( *j );
0228     trk2avg_[i] = pos.wtAverage(trks).value();
0229   }
0230 
0231   // Now let's send off our tracks to the Kalman fitter to see what great things happen...
0232   if (nvtx2_ > 0) {
0233     vector<reco::TransientTrack> t_tks;
0234     for (int i=0; i<nvtx2_ && i<maxvtx_; i++) {
0235       t_tks.clear();
0236       for (reco::track_iterator j=vertexes[i].tracks_begin(); j!=vertexes[i].tracks_end(); ++j) {
0237     t_tks.push_back( reco::TransientTrack(**j,field.product()) );
0238       }
0239       KalmanVertexFitter kvf;
0240       //      TransientVertex tv = kvf.vertex(t_tks,GlobalPoint(vertexes[i].x(),vertexes[i].y(),vertexes[i].z()));
0241       TransientVertex tv = kvf.vertex(t_tks);
0242       if (verbose_>0) std::cout << "Kalman Position: " << reco::Vertex::Point(tv.position()) << std::endl;
0243       vxkal_[i] = tv.position().x();
0244       vykal_[i] = tv.position().y();
0245       vzkal_[i] = tv.position().z();
0246     }
0247   }
0248     
0249   
0250 
0251   // Finally, fill the tree with the above values
0252   t_->Fill();
0253 }
0254 
0255 void PixelVertexTest::endJob() {
0256   if (t_) t_->Print();
0257   if (f_) {
0258     f_->Print();
0259     f_->Write();
0260   }
0261 }
0262 
0263 DEFINE_FWK_MODULE(PixelVertexTest);