File indexing completed on 2023-03-17 11:23:17
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
0155 DEFINE_FWK_MODULE(CVRAnalysis);
