TrackRefitter/test/TrajectoryReader.cc

0001 /** \class TrajectoryReader
0002  *  No description available.
0003  *
0004  *  \author R. Bellan - INFN Torino <riccardo.bellan@cern.ch>
0005  */
0006 // Base Class Headers
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   /// Constructor
0039   TrajectoryReader(const edm::ParameterSet &pset);
0040
0041   /// Destructor
0042   virtual ~TrajectoryReader();
0043
0044   void analyze(const edm::Event &event, const edm::EventSetup &eventSetup);
0045
0046   // Operations
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 /// Constructor
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 /// Destructor
0085 TrajectoryReader::~TrajectoryReader() = default;
0086
0087 // Operations
0088 void TrajectoryReader::beginJob() {
0089   // Create the root file
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   // Write the histos to file
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   // Global Tracking Geometry
0150   const GlobalTrackingGeometry *trackingGeometry = &eventSetup.getData(theGlobGeomToken);
0151
0152   // Magfield Field
0153   const MagneticField *magField = &eventSetup.getData(theMFToken);
0154
0155   const std::string metname = "Reco|TrackingTools|TrajectoryReader";
0156
0157   // Get the RecTrack collection from the event
0158   const reco::TrackCollection tracks = event.get(theTrackToken);
0159
0160   if (tracks.empty())
0161     return;
0162
0163   // Get the Trajectory collection from the event
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   //for (TrajTrackAssociationCollection::const_iterator it = assoMap.begin(); it != assoMap.end(); ++it) {
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     // Check the difference in Pt
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);