File indexing completed on 2024-04-06 12:31:39
0001
0002
0003
0004
0005
0006
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "DataFormats/DetId/interface/DetId.h"
0009 #include "DataFormats/TrackReco/interface/Track.h"
0010 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0011 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/Utilities/interface/InputTag.h"
0018 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0019 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0020 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0021 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0022 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0023 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0025 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0026 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0027
0028 #include "TFile.h"
0029 #include "TH1F.h"
0030
0031 #include <vector>
0032
0033 class TrajectoryReader : public edm::one::EDAnalyzer<> {
0034 public:
0035 typedef std::vector<Trajectory> Trajectories;
0036
0037 public:
0038
0039 TrajectoryReader(const edm::ParameterSet &pset);
0040
0041
0042 virtual ~TrajectoryReader();
0043
0044 void analyze(const edm::Event &event, const edm::EventSetup &eventSetup);
0045
0046
0047 void beginJob();
0048 void endJob();
0049
0050 protected:
0051 void printTrajectoryRecHits(const Trajectory &, edm::ESHandle<GlobalTrackingGeometry>) const;
0052 void printTrackRecHits(const reco::Track &, edm::ESHandle<GlobalTrackingGeometry>) const;
0053
0054 private:
0055 const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> theGlobGeomToken;
0056 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> theMFToken;
0057
0058 edm::EDGetTokenT<Trajectories> theTrajToken;
0059 edm::EDGetTokenT<reco::TrackCollection> theTrackToken;
0060 edm::EDGetTokenT<TrajTrackAssociationCollection> theAssocToken;
0061
0062 edm::InputTag theInputLabel;
0063 TFile *theFile;
0064 std::string theRootFileName;
0065
0066 TH1F *hDPtIn;
0067 TH1F *hDPtOut;
0068 TH1F *hSuccess;
0069 TH1F *hNHitLost;
0070 TH1F *hFractionHitLost;
0071 };
0072
0073
0074 TrajectoryReader::TrajectoryReader(const edm::ParameterSet &pset)
0075 : theGlobGeomToken(esConsumes()), theMFToken(esConsumes()) {
0076 theInputLabel = pset.getParameter<edm::InputTag>("InputLabel");
0077
0078 theTrajToken = consumes<Trajectories>(theInputLabel), theTrackToken = consumes<reco::TrackCollection>(theInputLabel),
0079 theAssocToken = consumes<TrajTrackAssociationCollection>(theInputLabel);
0080
0081 theRootFileName = pset.getUntrackedParameter<std::string>("rootFileName");
0082 }
0083
0084
0085 TrajectoryReader::~TrajectoryReader() = default;
0086
0087
0088 void TrajectoryReader::beginJob() {
0089
0090 theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0091 theFile->cd();
0092
0093 hDPtIn = new TH1F("DeltaPtIn", "P_{t}^{Track}-P_{t}^{Traj} inner state", 10000, -20, 20);
0094 hDPtOut = new TH1F("DeltaPtOut", "P_{t}^{Track}-P_{t}^{Traj} outer state", 10000, -20, 20);
0095
0096 hNHitLost = new TH1F("NHitLost", "Number of lost hits", 100, 0, 100);
0097 hFractionHitLost = new TH1F("FractionHitLost", "Fraction of lost hits", 100, 0, 100);
0098 hSuccess = new TH1F("Success", "Number of Success", 2, 0, 2);
0099 }
0100
0101 void TrajectoryReader::endJob() {
0102 theFile->cd();
0103
0104
0105 hDPtIn->Write();
0106 hDPtOut->Write();
0107 hNHitLost->Write();
0108 hFractionHitLost->Write();
0109 hSuccess->Write();
0110 theFile->Close();
0111 }
0112
0113 void TrajectoryReader::printTrajectoryRecHits(const Trajectory &trajectory,
0114 edm::ESHandle<GlobalTrackingGeometry> trackingGeometry) const {
0115 const std::string metname = "Reco|TrackingTools|TrajectoryReader";
0116
0117 TransientTrackingRecHit::ConstRecHitContainer rechits = trajectory.recHits();
0118 LogDebug(metname) << "Size of the RecHit container: " << rechits.size();
0119
0120 int i = 0;
0121 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator recHit = rechits.begin(); recHit != rechits.end();
0122 ++recHit)
0123 if ((*recHit)->isValid()) {
0124 const GeomDet *geomDet = trackingGeometry->idToDet((*recHit)->geographicalId());
0125 double r = geomDet->surface().position().perp();
0126 double z = geomDet->toGlobal((*recHit)->localPosition()).z();
0127 LogTrace(metname) << i++ << " r: " << r << " z: " << z << " " << geomDet->toGlobal((*recHit)->localPosition())
0128 << std::endl;
0129 }
0130 }
0131
0132 void TrajectoryReader::printTrackRecHits(const reco::Track &track,
0133 edm::ESHandle<GlobalTrackingGeometry> trackingGeometry) const {
0134 const std::string metname = "Reco|TrackingTools|TrajectoryReader";
0135
0136 LogTrace(metname) << "Valid RecHits: " << track.found() << " invalid RecHits: " << track.lost();
0137
0138 int i = 0;
0139 for (trackingRecHit_iterator recHit = track.recHitsBegin(); recHit != track.recHitsEnd(); ++recHit)
0140 if ((*recHit)->isValid()) {
0141 const GeomDet *geomDet = trackingGeometry->idToDet((*recHit)->geographicalId());
0142 double r = geomDet->surface().position().perp();
0143 double z = geomDet->surface().position().z();
0144 LogTrace(metname) << i++ << " GeomDet position r: " << r << " z: " << z;
0145 }
0146 }
0147
0148 void TrajectoryReader::analyze(const edm::Event &event, const edm::EventSetup &eventSetup) {
0149
0150 const GlobalTrackingGeometry *trackingGeometry = &eventSetup.getData(theGlobGeomToken);
0151
0152
0153 const MagneticField *magField = &eventSetup.getData(theMFToken);
0154
0155 const std::string metname = "Reco|TrackingTools|TrajectoryReader";
0156
0157
0158 const reco::TrackCollection tracks = event.get(theTrackToken);
0159
0160 if (tracks.empty())
0161 return;
0162
0163
0164 const Trajectories trajectories = event.get(theTrajToken);
0165
0166 LogTrace(metname) << "looking at: " << theInputLabel;
0167
0168 LogTrace(metname) << "All trajectories";
0169 for (const auto &trajectory : trajectories) {
0170 printTrajectoryRecHits(trajectory, trackingGeometry);
0171 }
0172
0173 LogTrace(metname) << "All tracks";
0174 for (const auto &tr : tracks) {
0175 printTrackRecHits(tr, trackingGeometry);
0176 }
0177
0178 const TrajTrackAssociationCollection assoMap = event.get(theAssocToken);
0179
0180 LogTrace(metname) << "Association";
0181
0182 for (const auto &it : assoMap) {
0183 const edm::Ref<std::vector<Trajectory> > traj = it.key;
0184 const reco::TrackRef tk = it.val;
0185
0186 printTrackRecHits(*tk, trackingGeometry);
0187 printTrajectoryRecHits(*traj, trackingGeometry);
0188
0189
0190 reco::TransientTrack track(tk, &*magField, trackingGeometry);
0191
0192 hDPtIn->Fill(track.innermostMeasurementState().globalMomentum().perp() -
0193 traj->lastMeasurement().updatedState().globalMomentum().perp());
0194 hDPtOut->Fill(track.outermostMeasurementState().globalMomentum().perp() -
0195 traj->firstMeasurement().updatedState().globalMomentum().perp());
0196
0197 int diff = track.recHitsSize() - traj->recHits().size();
0198 LogTrace(metname) << "Difference: " << diff;
0199 hNHitLost->Fill(diff);
0200 hFractionHitLost->Fill(double(diff) / track.recHitsSize());
0201 }
0202
0203 int traj_size = trajectories.size();
0204 int track_size = tracks.size();
0205
0206 if (traj_size != track_size) {
0207 LogTrace(metname) << "Mismatch between the # of Tracks (" << track_size << ") and the # of Trajectories! ("
0208 << traj_size << ") in " << event.id();
0209 hSuccess->Fill(0);
0210 } else
0211 hSuccess->Fill(1);
0212 }
0213
0214 #include "FWCore/PluginManager/interface/ModuleDef.h"
0215 #include "FWCore/Framework/interface/MakerMacros.h"
0216
0217 DEFINE_FWK_MODULE(TrajectoryReader);