DebugTools/plugins/TrackAlgoCompareUtil.cc

0001 #include "RecoTracker/DebugTools/plugins/TrackAlgoCompareUtil.h"
0002
0003 using namespace std;
0004 using namespace edm;
0005
0006 // constructors and destructor
0007 TrackAlgoCompareUtil::TrackAlgoCompareUtil(const edm::ParameterSet& iConfig)
0008     : trackLabel_algoA(consumes<View<reco::Track>>(iConfig.getParameter<edm::InputTag>("trackLabel_algoA"))),
0009       trackLabel_algoB(consumes<View<reco::Track>>(iConfig.getParameter<edm::InputTag>("trackLabel_algoB"))),
0010       trackingParticleLabel_fakes(
0011           consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("trackingParticleLabel_fakes"))),
0012       trackingParticleLabel_effic(
0013           consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("trackingParticleLabel_effic"))),
0014       beamSpotLabel(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpotLabel"))),
0015       magField(esConsumes()),
0016       UseAssociators(iConfig.getParameter<bool>("UseAssociators")),
0017       UseVertex(iConfig.getParameter<bool>("UseVertex")) {
0018   //now do what ever other initialization is needed
0019   if (UseVertex) {
0020     vertexLabel_algoA = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexLabel_algoA"));
0021     vertexLabel_algoB = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexLabel_algoB"));
0022   }
0023
0024   if (UseAssociators) {
0025     assocLabel_algoA = consumes<reco::TrackToTrackingParticleAssociator>(
0026         iConfig.getUntrackedParameter<std::string>("assocLabel_algoA", "trackAssociatorByHits"));
0027     assocLabel_algoB = consumes<reco::TrackToTrackingParticleAssociator>(
0028         iConfig.getUntrackedParameter<std::string>("assocLabel_algoB", "trackAssociatorByHits"));
0029   } else {
0030     edm::InputTag algoA = iConfig.getParameter<edm::InputTag>("associatormap_algoA");
0031     edm::InputTag algoB = iConfig.getParameter<edm::InputTag>("associatormap_algoB");
0032
0033     associatormap_algoA_recoToSim = consumes<reco::RecoToSimCollection>(algoA);
0034     associatormap_algoB_recoToSim = consumes<reco::RecoToSimCollection>(algoB);
0035     associatormap_algoA_simToReco = consumes<reco::SimToRecoCollection>(algoA);
0036     associatormap_algoB_simToReco = consumes<reco::SimToRecoCollection>(algoB);
0037   }
0038
0039   produces<RecoTracktoTPCollection>("AlgoA");
0040   produces<RecoTracktoTPCollection>("AlgoB");
0041   produces<TPtoRecoTrackCollection>("TP");
0042 }
0043
0044 TrackAlgoCompareUtil::~TrackAlgoCompareUtil() {}
0045
0046 // ------------ method called to produce the data  ------------
0047 void TrackAlgoCompareUtil::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0048   // create output collection instance
0049   auto outputAlgoA = std::make_unique<RecoTracktoTPCollection>();
0050   auto outputAlgoB = std::make_unique<RecoTracktoTPCollection>();
0051   auto outputTP = std::make_unique<TPtoRecoTrackCollection>();
0052
0053   // Get Inputs
0054   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0055   iEvent.getByToken(beamSpotLabel, recoBeamSpotHandle);
0056   reco::BeamSpot beamSpot = *recoBeamSpotHandle;
0057
0058   edm::Handle<View<reco::Track>> trackCollAlgoA;
0059   iEvent.getByToken(trackLabel_algoA, trackCollAlgoA);
0060
0061   edm::Handle<View<reco::Track>> trackCollAlgoB;
0062   iEvent.getByToken(trackLabel_algoB, trackCollAlgoB);
0063
0064   edm::Handle<TrackingParticleCollection> trackingParticleCollFakes;
0065   iEvent.getByToken(trackingParticleLabel_fakes, trackingParticleCollFakes);
0066
0067   edm::Handle<TrackingParticleCollection> trackingParticleCollEffic;
0068   iEvent.getByToken(trackingParticleLabel_effic, trackingParticleCollEffic);
0069
0070   edm::Handle<reco::VertexCollection> vertexCollAlgoA;
0071   edm::Handle<reco::VertexCollection> vertexCollAlgoB;
0072   if (UseVertex) {
0073     iEvent.getByToken(vertexLabel_algoA, vertexCollAlgoA);
0074     iEvent.getByToken(vertexLabel_algoB, vertexCollAlgoB);
0075   }
0076
0077   // call the associator functions:
0078   reco::RecoToSimCollection recSimColl_AlgoA;
0079   reco::RecoToSimCollection recSimColl_AlgoB;
0080
0081   reco::SimToRecoCollection simRecColl_AlgoA;
0082   reco::SimToRecoCollection simRecColl_AlgoB;
0083
0084   if (UseAssociators) {
0085     edm::Handle<reco::TrackToTrackingParticleAssociator> theAssociator_algoA;
0086     iEvent.getByToken(assocLabel_algoA, theAssociator_algoA);
0087
0088     edm::Handle<reco::TrackToTrackingParticleAssociator> theAssociator_algoB;
0089     iEvent.getByToken(assocLabel_algoB, theAssociator_algoB);
0090
0091     recSimColl_AlgoA = theAssociator_algoA->associateRecoToSim(trackCollAlgoA, trackingParticleCollFakes);
0092     recSimColl_AlgoB = theAssociator_algoB->associateRecoToSim(trackCollAlgoB, trackingParticleCollFakes);
0093
0094     simRecColl_AlgoA = theAssociator_algoA->associateSimToReco(trackCollAlgoA, trackingParticleCollEffic);
0095     simRecColl_AlgoB = theAssociator_algoB->associateSimToReco(trackCollAlgoB, trackingParticleCollEffic);
0096   } else {
0097     Handle<reco::RecoToSimCollection> recotosimCollectionH_AlgoA;
0098     iEvent.getByToken(associatormap_algoA_recoToSim, recotosimCollectionH_AlgoA);
0099     recSimColl_AlgoA = *(recotosimCollectionH_AlgoA.product());
0100
0101     Handle<reco::RecoToSimCollection> recotosimCollectionH_AlgoB;
0102     iEvent.getByToken(associatormap_algoB_recoToSim, recotosimCollectionH_AlgoB);
0103     recSimColl_AlgoB = *(recotosimCollectionH_AlgoB.product());
0104
0105     Handle<reco::SimToRecoCollection> simtorecoCollectionH_AlgoA;
0106     iEvent.getByToken(associatormap_algoA_simToReco, simtorecoCollectionH_AlgoA);
0107     simRecColl_AlgoA = *(simtorecoCollectionH_AlgoA.product());
0108
0109     Handle<reco::SimToRecoCollection> simtorecoCollectionH_AlgoB;
0110     iEvent.getByToken(associatormap_algoB_simToReco, simtorecoCollectionH_AlgoB);
0111     simRecColl_AlgoB = *(simtorecoCollectionH_AlgoB.product());
0112   }
0113
0114   // define the vector of references to trackingParticleColl associated with a given reco::Track
0115   std::vector<std::pair<TrackingParticleRef, double>> associatedTrackingParticles;
0116
0117   // define the vector of references to trackColl associated with a given TrackingParticle
0118   std::vector<std::pair<reco::TrackBaseRef, double>> associatedRecoTracks;
0119
0120   // Get the magnetic field data from the event (used to calculate the point of closest TrackingParticle)
0121   const MagneticField* magneticField = &iSetup.getData(magField);
0122
0123   // fill collection algoA
0124   for (View<reco::Track>::size_type i = 0; i < trackCollAlgoA->size(); ++i) {
0125     // get recoTrack algo A
0126     reco::TrackBaseRef recoTrack(trackCollAlgoA, i);
0127     RecoTracktoTP recoTracktoTP;
0128     recoTracktoTP.SetRecoTrack(recoTrack);
0129     recoTracktoTP.SetBeamSpot(beamSpot.position());
0130
0131     // get the associated trackingParticle
0132     if (recSimColl_AlgoA.find(recoTrack) != recSimColl_AlgoA.end()) {
0133       associatedTrackingParticles = recSimColl_AlgoA[recoTrack];
0134       recoTracktoTP.SetTrackingParticle(associatedTrackingParticles.begin()->first);
0135       recoTracktoTP.SetShared(associatedTrackingParticles.begin()->second);
0136       SetTrackingParticleD0Dz(associatedTrackingParticles.begin()->first, beamSpot, magneticField, recoTracktoTP);
0137     } else {
0138       recoTracktoTP.SetTrackingParticle(TrackingParticleRef());
0139       recoTracktoTP.SetShared(-1.0);
0140     }
0141
0142     // get the reco primary vertex info
0143     if (UseVertex && !vertexCollAlgoA->empty()) {
0144       recoTracktoTP.SetRecoVertex(reco::VertexRef(vertexCollAlgoA, 0));
0145     } else {
0146       recoTracktoTP.SetRecoVertex(reco::VertexRef());
0147     }
0148
0149     outputAlgoA->push_back(recoTracktoTP);
0150   }
0151
0152   // fill collection algoB
0153   for (reco::TrackCollection::size_type i = 0; i < trackCollAlgoB->size(); ++i) {
0154     // get recoTrack algo B
0155     reco::TrackBaseRef recoTrack(trackCollAlgoB, i);
0156     RecoTracktoTP recoTracktoTP;
0157     recoTracktoTP.SetRecoTrack(recoTrack);
0158     recoTracktoTP.SetBeamSpot(beamSpot.position());
0159
0160     // get the associated trackingParticle
0161     if (recSimColl_AlgoB.find(recoTrack) != recSimColl_AlgoB.end()) {
0162       associatedTrackingParticles = recSimColl_AlgoB[recoTrack];
0163       recoTracktoTP.SetTrackingParticle(associatedTrackingParticles.begin()->first);
0164       recoTracktoTP.SetShared(associatedTrackingParticles.begin()->second);
0165       SetTrackingParticleD0Dz(associatedTrackingParticles.begin()->first, beamSpot, magneticField, recoTracktoTP);
0166     } else {
0167       recoTracktoTP.SetTrackingParticle(TrackingParticleRef());
0168       recoTracktoTP.SetShared(-1.0);
0169     }
0170
0171     // get the reco primary vertex info
0172     if (UseVertex && !vertexCollAlgoB->empty()) {
0173       recoTracktoTP.SetRecoVertex(reco::VertexRef(vertexCollAlgoB, 0));
0174     } else {
0175       recoTracktoTP.SetRecoVertex(reco::VertexRef());
0176     }
0177
0178     outputAlgoB->push_back(recoTracktoTP);
0179   }
0180
0181   for (TrackingParticleCollection::size_type i = 0; i < trackingParticleCollEffic->size(); ++i) {
0182     // initialize the trackingParticle (sim) info
0183     TrackingParticleRef tparticle(trackingParticleCollEffic, i);
0184     TPtoRecoTrack tptoRecoTrack;
0185     tptoRecoTrack.SetBeamSpot(beamSpot.position());
0186     tptoRecoTrack.SetTrackingParticle(tparticle);
0187     SetTrackingParticleD0Dz(tparticle, beamSpot, magneticField, tptoRecoTrack);
0188
0189     // get the assocated recoTrack algoA
0190     if (simRecColl_AlgoA.find(tparticle) != simRecColl_AlgoA.end()) {
0191       associatedRecoTracks = simRecColl_AlgoA[tparticle];
0192       tptoRecoTrack.SetRecoTrack_AlgoA(associatedRecoTracks.begin()->first);
0193       tptoRecoTrack.SetShared_AlgoA(associatedRecoTracks.begin()->second);
0194     } else {
0195       tptoRecoTrack.SetRecoTrack_AlgoA(reco::TrackBaseRef());
0196       tptoRecoTrack.SetShared_AlgoA(-1.0);
0197     }
0198
0199     // get the recoVertex algo A
0200     if (UseVertex && !vertexCollAlgoA->empty()) {
0201       tptoRecoTrack.SetRecoVertex_AlgoA(reco::VertexRef(vertexCollAlgoA, 0));
0202     } else {
0203       tptoRecoTrack.SetRecoVertex_AlgoA(reco::VertexRef());
0204     }
0205
0206     // get the assocated recoTrack algoB
0207     if (simRecColl_AlgoB.find(tparticle) != simRecColl_AlgoB.end()) {
0208       associatedRecoTracks = simRecColl_AlgoB[tparticle];
0209       tptoRecoTrack.SetRecoTrack_AlgoB(associatedRecoTracks.begin()->first);
0210       tptoRecoTrack.SetShared_AlgoB(associatedRecoTracks.begin()->second);
0211     } else {
0212       tptoRecoTrack.SetRecoTrack_AlgoB(reco::TrackBaseRef());
0213       tptoRecoTrack.SetShared_AlgoB(-1.0);
0214     }
0215     // get the recoVertex algo B
0216     if (UseVertex && !vertexCollAlgoB->empty()) {
0217       tptoRecoTrack.SetRecoVertex_AlgoB(reco::VertexRef(vertexCollAlgoB, 0));
0218     } else {
0219       tptoRecoTrack.SetRecoVertex_AlgoB(reco::VertexRef());
0220     }
0221
0222     outputTP->push_back(tptoRecoTrack);
0223   }
0224
0225   // put the collection in the event record
0226   iEvent.put(std::move(outputAlgoA), "AlgoA");
0227   iEvent.put(std::move(outputAlgoB), "AlgoB");
0228   iEvent.put(std::move(outputTP), "TP");
0229 }
0230
0231 // ------------ Producer Specific Meber Fucntions ----------------------------------------
0232 void TrackAlgoCompareUtil::SetTrackingParticleD0Dz(TrackingParticleRef tp,
0233                                                    const reco::BeamSpot& bs,
0234                                                    const MagneticField* bf,
0235                                                    TPtoRecoTrack& TPRT) const {
0236   GlobalPoint trackingParticleVertex(tp->vertex().x(), tp->vertex().y(), tp->vertex().z());
0237   GlobalVector trackingParticleP3(
0238       tp->g4Track_begin()->momentum().x(), tp->g4Track_begin()->momentum().y(), tp->g4Track_begin()->momentum().z());
0239   TrackCharge trackingParticleCharge(tp->charge());
0240
0241   FreeTrajectoryState ftsAtProduction(trackingParticleVertex, trackingParticleP3, trackingParticleCharge, bf);
0242   TSCBLBuilderNoMaterial tscblBuilder;
0243   TrajectoryStateClosestToBeamLine tsAtClosestApproach =
0244       tscblBuilder(ftsAtProduction, bs);  //as in TrackProducerAlgorithm
0245
0246   if (tsAtClosestApproach.isValid()) {
0247     GlobalPoint v1 = tsAtClosestApproach.trackStateAtPCA().position();
0248     GlobalVector p = tsAtClosestApproach.trackStateAtPCA().momentum();
0249
0250     TPRT.SetTrackingParticleMomentumPCA(p);
0251     TPRT.SetTrackingParticlePCA(v1);
0252   } else {
0253     TPRT.SetTrackingParticleMomentumPCA(GlobalVector(-9999.0, -9999.0, -9999.0));
0254     TPRT.SetTrackingParticlePCA(GlobalPoint(-9999.0, -9999.0, -9999.0));
0255   }
0256 }
0257
0258 void TrackAlgoCompareUtil::SetTrackingParticleD0Dz(TrackingParticleRef tp,
0259                                                    const reco::BeamSpot& bs,
0260                                                    const MagneticField* bf,
0261                                                    RecoTracktoTP& RTTP) const {
0262   GlobalPoint trackingParticleVertex(tp->vertex().x(), tp->vertex().y(), tp->vertex().z());
0263   GlobalVector trackingParticleP3(
0264       tp->g4Track_begin()->momentum().x(), tp->g4Track_begin()->momentum().y(), tp->g4Track_begin()->momentum().z());
0265   TrackCharge trackingParticleCharge(tp->charge());
0266
0267   FreeTrajectoryState ftsAtProduction(trackingParticleVertex, trackingParticleP3, trackingParticleCharge, bf);
0268   TSCBLBuilderNoMaterial tscblBuilder;
0269   TrajectoryStateClosestToBeamLine tsAtClosestApproach =
0270       tscblBuilder(ftsAtProduction, bs);  //as in TrackProducerAlgorithm
0271
0272   if (tsAtClosestApproach.isValid()) {
0273     GlobalPoint v1 = tsAtClosestApproach.trackStateAtPCA().position();
0274     GlobalVector p = tsAtClosestApproach.trackStateAtPCA().momentum();
0275
0276     RTTP.SetTrackingParticleMomentumPCA(p);
0277     RTTP.SetTrackingParticlePCA(v1);
0278   } else {
0279     RTTP.SetTrackingParticleMomentumPCA(GlobalVector(-9999.0, -9999.0, -9999.0));
0280     RTTP.SetTrackingParticlePCA(GlobalPoint(-9999.0, -9999.0, -9999.0));
0281   }
0282 }