Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:28:08

0001 #include <vector>
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Framework/interface/ESHandle.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "MagneticField/Engine/interface/MagneticField.h"
0007 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0008 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0009 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0010 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0011 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0012 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0013 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0014 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0015 #include "RecoVertex/ConfigurableVertexReco/test/CVRAnalysis.h"
0016 #include <iostream>
0017 
0018 using namespace std;
0019 using namespace reco;
0020 using namespace edm;
0021 
0022 
0023 namespace {
0024   void printTSOS ( const TrajectoryStateOnSurface & tsos )
0025   {
0026     cout << tsos.globalPosition() << " , "
0027          << tsos.globalMomentum() << endl;
0028   }
0029 
0030   void printVertex ( const TransientVertex & vtx )
0031   {
0032     cout << " `- pos=(" << vtx.position().x() << ", "
0033          << vtx.position().y() << ", " << vtx.position().z()
0034          << ") chi2=" << vtx.totalChiSquared()
0035          << " ndf=" << vtx.degreesOfFreedom() << " hr="
0036          << vtx.hasRefittedTracks() << endl;
0037     if ( vtx.originalTracks().size() && vtx.hasRefittedTracks() )
0038     {
0039       cout << "    `- 1st trk: ";
0040       reco::TransientTrack t = vtx.originalTracks()[0];
0041       TrajectoryStateOnSurface tsos = t.impactPointState();
0042       printTSOS ( tsos );
0043       if ( vtx.refittedTracks().size() )
0044       {
0045         cout << "     `- 1st refttd: ";
0046         reco::TransientTrack t2 = vtx.refittedTracks()[0];
0047         printTSOS ( t2.impactPointState() );
0048       }
0049     }
0050   }
0051 
0052   void printVertices ( const vector < TransientVertex > & vtces )
0053   {
0054     cout << "[CVRAnalysis] " << vtces.size() << " vertices." << endl;
0055     for ( vector< TransientVertex >::const_iterator i=vtces.begin();
0056           i!=vtces.end() ; ++i )
0057     {
0058       printVertex ( *i );
0059       cout << endl;
0060     }
0061   }
0062 
0063   void discussBeamSpot ( const reco::BeamSpot & bs )
0064   {
0065     cout << "[CVRAnalysis] beamspot at " << bs.position() << endl;
0066     reco::BeamSpot::Covariance3DMatrix cov = bs.rotatedCovariance3D();
0067     cout << "[CVRAnalysis] cov=" <<  cov << endl;
0068   }
0069 }
0070 
0071 CVRAnalysis::CVRAnalysis(const edm::ParameterSet& iconfig) :
0072   trackcoll_( iconfig.getParameter<string>("trackcoll") ),
0073   vertexcoll_( iconfig.getParameter<string>("vertexcoll") ),
0074   beamspot_( iconfig.getParameter<string>("beamspot") ),
0075   trackingtruth_ ( iconfig.getParameter< edm::InputTag >("truth") ),
0076   associator_ ( iconfig.getParameter<string>("associator") ),
0077   histo_ ( VertexHisto ( "vertices.root", "tracks.root" ) ),
0078   bhisto_ ( VertexHisto ( "vertices-b.root", "tracks-b.root" ) )
0079 {
0080   edm::ParameterSet vtxconfig = iconfig.getParameter<edm::ParameterSet>("vertexreco");
0081   vrec_ = new ConfigurableVertexReconstructor ( vtxconfig );
0082   cout << "[CVRAnalysis] vtxconfig=" << vtxconfig << endl;
0083 }
0084 
0085 CVRAnalysis::~CVRAnalysis() {
0086   if ( vrec_ ) delete vrec_;
0087 }
0088 
0089 void CVRAnalysis::discussPrimary( const edm::Event& iEvent ) const
0090 {
0091   edm::Handle<reco::VertexCollection> retColl;
0092   iEvent.getByLabel( vertexcoll_, retColl);
0093   if ( retColl->size() )
0094   {
0095     const reco::Vertex & vtx = *(retColl->begin());
0096     cout << "[CVRAnalysis] persistent primary: " << vtx.x() << ", " << vtx.y()
0097          << ", " << vtx.z() << endl;
0098   }
0099 }
0100 
0101 void CVRAnalysis::analyze( const edm::Event & iEvent,
0102                        const edm::EventSetup & iSetup )
0103 {
0104   int evt=iEvent.id().event();
0105   cout << "[CVRAnalysis] next event: " << evt << endl;
0106   edm::ESHandle<MagneticField> magneticField;
0107   iSetup.get<IdealMagneticFieldRecord>().get(magneticField);
0108   edm::ESHandle<TransientTrackBuilder> builder;
0109   iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",builder );
0110 
0111   edm::Handle< edm::View < reco::Track > > tks;
0112   iEvent.getByLabel( trackcoll_, tks );
0113   discussPrimary( iEvent );
0114 
0115   edm::Handle<TrackingVertexCollection> TVCollectionH;
0116   iEvent.getByLabel( trackingtruth_, TVCollectionH);
0117 
0118   edm::Handle<TrackingParticleCollection>  TPCollectionH;
0119   iEvent.getByLabel( trackingtruth_, TPCollectionH);
0120 
0121   edm::Handle<reco::TrackToTrackingParticleAssociator> byHitsAssociator;
0122   iEvent.getByLabel(associator_, byHitsAssociator);
0123 
0124   reco::RecoToSimCollection p =
0125   byHitsAssociator->associateRecoToSim ( tks, TPCollectionH );
0126 
0127   edm::Handle<reco::BeamSpot > bs;
0128   iEvent.getByLabel ( beamspot_, bs );
0129   discussBeamSpot ( *bs );
0130 
0131   vector<reco::TransientTrack> ttks;
0132   ttks = builder->build(tks);
0133   cout << "[CVRAnalysis] got " << ttks.size() << " tracks." << endl;
0134 
0135   cout << "[CVRAnalysis] fit w/o beamspot constraint" << endl;
0136   vector < TransientVertex > vtces = vrec_->vertices ( ttks );
0137   printVertices ( vtces );
0138 
0139   if ( vtces.size() && TVCollectionH->size() )
0140   {
0141     histo_.analyse ( *(TVCollectionH->begin()), vtces[0], "Primaries" );
0142     histo_.saveTracks ( vtces[0], p, "VtxTk" );
0143   }
0144 
0145   cout << "[CVRAnalysis] fit w beamspot constraint" << endl;
0146   vector < TransientVertex > bvtces = vrec_->vertices ( ttks, *bs );
0147   printVertices ( bvtces );
0148   if ( bvtces.size() && TVCollectionH->size() )
0149   {
0150     bhisto_.analyse ( *(TVCollectionH->begin()), bvtces[0], "Primaries" );
0151   }
0152 }
0153 
0154 //define this as a plug-in
0155 DEFINE_FWK_MODULE(CVRAnalysis);