TrackingMonitor/src/TrackingRecoMaterialAnalyser.cc

0001 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0002 #include "DataFormats/TrackReco/interface/Track.h"
0003 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0004 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0005 #include "DataFormats/DetId/interface/DetId.h"
0006 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0007 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0008 #include "DataFormats/VertexReco/interface/Vertex.h"
0009 #include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h"
0010
0011 #include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"
0012 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0013 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0014 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0015 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0016
0017 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0018 #include "DQMServices/Core/interface/DQMStore.h"
0019 #include "FWCore/Framework/interface/ConsumesCollector.h"
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023
0024 #include <cassert>
0025 #include <unordered_map>
0026 #include <string>
0027
0028 // Values are not ordered randomly, but the order is taken from
0029 // http://cmslxr.fnal.gov/dxr/CMSSW/source/Geometry/CommonDetUnit/interface/GeomDetEnumerators.h#15
0030 static const std::vector<std::string> sDETS{"", "PXB", "PXF", "TIB", "TID", "TOB", "TEC"};
0031
0032 class TrackingRecoMaterialAnalyser : public DQMEDAnalyzer {
0033 public:
0034   explicit TrackingRecoMaterialAnalyser(const edm::ParameterSet &);
0035   void bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) override;
0036   void analyze(const edm::Event &, const edm::EventSetup &) override;
0037   ~TrackingRecoMaterialAnalyser() override;
0038
0039 private:
0040   inline bool isDoubleSided(const TrackerTopology *tTopo, DetId id) { return (tTopo->glued(id)); }
0041   TrackTransformer refitter_;
0042   const edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0043   const edm::EDGetTokenT<reco::BeamSpot> beamspotToken_;
0044   const edm::EDGetTokenT<reco::VertexCollection> verticesToken_;
0045   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeometryTokenRun_;
0046   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0047   bool usePV_;
0048   std::string folder_;
0049   std::unordered_map<std::string, MonitorElement *> histosOriEta_;
0050   std::unordered_map<std::string, MonitorElement *> histosEta_;
0051   MonitorElement *histo_RZ_;
0052   MonitorElement *histo_RZ_Ori_;
0053   MonitorElement *deltaPt_in_out_2d_;
0054   MonitorElement *deltaP_in_out_vs_eta_;
0055   MonitorElement *deltaP_in_out_vs_z_;
0056   MonitorElement *deltaP_in_out_vs_eta_2d_;
0057   MonitorElement *deltaP_in_out_vs_eta_vs_phi_2d_;
0058   MonitorElement *deltaP_in_out_vs_z_2d_;
0059   MonitorElement *deltaPt_in_out_vs_eta_;
0060   MonitorElement *deltaPt_in_out_vs_z_;
0061   MonitorElement *deltaPl_in_out_vs_eta_;
0062   MonitorElement *deltaPl_in_out_vs_z_;
0063   MonitorElement *P_vs_eta_2d_;
0064 };
0065
0066 //-------------------------------------------------------------------------
0067 TrackingRecoMaterialAnalyser::TrackingRecoMaterialAnalyser(const edm::ParameterSet &iPSet)
0068     : refitter_(iPSet, consumesCollector()),
0069       tracksToken_(consumes<reco::TrackCollection>(iPSet.getParameter<edm::InputTag>("tracks"))),
0070       beamspotToken_(consumes<reco::BeamSpot>(iPSet.getParameter<edm::InputTag>("beamspot"))),
0071       verticesToken_(mayConsume<reco::VertexCollection>(iPSet.getParameter<edm::InputTag>("vertices"))),
0072       trackerGeometryTokenRun_(esConsumes<edm::Transition::BeginRun>()),
0073       tTopoToken_(esConsumes()),
0074       usePV_(iPSet.getParameter<bool>("usePV")),
0075       folder_(iPSet.getParameter<std::string>("folder")),
0076       histo_RZ_(nullptr),
0077       histo_RZ_Ori_(nullptr),
0078       deltaPt_in_out_2d_(nullptr),
0079       deltaP_in_out_vs_eta_(nullptr),
0080       deltaP_in_out_vs_z_(nullptr),
0081       deltaP_in_out_vs_eta_2d_(nullptr),
0082       deltaP_in_out_vs_eta_vs_phi_2d_(nullptr),
0083       deltaP_in_out_vs_z_2d_(nullptr),
0084       deltaPt_in_out_vs_eta_(nullptr),
0085       deltaPt_in_out_vs_z_(nullptr),
0086       deltaPl_in_out_vs_eta_(nullptr),
0087       deltaPl_in_out_vs_z_(nullptr),
0088       P_vs_eta_2d_(nullptr) {}
0089
0090 //-------------------------------------------------------------------------
0091 TrackingRecoMaterialAnalyser::~TrackingRecoMaterialAnalyser(void) {}
0092
0093 //-------------------------------------------------------------------------
0094 void TrackingRecoMaterialAnalyser::bookHistograms(DQMStore::IBooker &ibook,
0095                                                   edm::Run const &,
0096                                                   edm::EventSetup const &setup) {
0097   using namespace std;
0098   const TrackerGeometry &trackerGeometry = setup.getData(trackerGeometryTokenRun_);
0099
0100   ibook.setCurrentFolder(folder_);
0101
0102   // Histogram to store the radiation length map, in the R-Z plane,
0103   // gathering numbers directly from the trackerRecoMaterial.xml
0104   // file. The numbers are not corrected for the track angle.
0105   histo_RZ_Ori_ = ibook.bookProfile2D("OriRadLen", "Original_RadLen", 600, -300., 300, 120, 0., 120., 0., 1.);
0106
0107   // Histogram to store the radiation length map, as before, but
0108   // correcting the numbers with the track angle. This represents the
0109   // real material seen by the track.
0110   histo_RZ_ = ibook.bookProfile2D("RadLen", "RadLen", 600, -300., 300, 120, 0., 120., 0., 1.);
0111
0112   // Histogram to show the deltaP Out-In in the eta-phi plane.
0113   deltaP_in_out_vs_eta_vs_phi_2d_ = ibook.bookProfile2D(
0114       "DeltaP_in_out_vs_eta_vs_phi_2d", "DeltaP_in_out_vs_eta_vs_phi_2d", 100, -3.0, 3.0, 100, -3.15, 3.15, 0., 100.);
0115
0116   // Histogram to show the deltaP Out-In vs eta.
0117   deltaP_in_out_vs_eta_2d_ =
0118       ibook.book2D("DeltaP_in_out_vs_eta_2d", "DeltaP_in_out_vs_eta_2d", 100, -3.0, 3.0, 100, 0., 1);
0119
0120   // Histogram to show the deltaP Out-In vs Z. The Z coordinate is the
0121   // one computed at the outermost hit.
0122   deltaP_in_out_vs_z_2d_ = ibook.book2D("DeltaP_in_out_vs_z_2d", "DeltaP_in_out_vs_z_2d", 600, -300, 300, 200., -1, 1.);
0123
0124   // Histogram to show the deltaP Out-In vs eta. The eta is the one
0125   // computed at the innermost hit.
0126   deltaP_in_out_vs_eta_ =
0127       ibook.bookProfile("DeltaP_in_out_vs_eta", "DeltaP_in_out_vs_eta", 100, -3.0, 3.0, -100., 100.);
0128   deltaP_in_out_vs_z_ = ibook.bookProfile("DeltaP_in_out_vs_z", "DeltaP_in_out_vs_z", 600, -300, 300, -100., 100.);
0129
0130   // Histogram to show the delta_Pt Out-In vs eta. The eta is the one
0131   // computed at the innermost hit.
0132   deltaPt_in_out_vs_eta_ =
0133       ibook.bookProfile("DeltaPt_in_out_vs_eta", "DeltaPt_in_out_vs_eta", 100, -3.0, 3.0, -100., 100.);
0134
0135   // Histogram to show the delta_Pt Out-In vs Z. The Z is the one
0136   // computed at the outermost hit.
0137   deltaPt_in_out_vs_z_ = ibook.bookProfile("DeltaPt_in_out_vs_z", "DeltaPt_in_out_vs_z", 600, -300, 300, -100., 100);
0138
0139   // Histogram to show the delta_Pl Out-In vs eta. The eta is the one
0140   // computed at the innermost hit.
0141   deltaPl_in_out_vs_eta_ =
0142       ibook.bookProfile("DeltaPz_in_out_vs_eta", "DeltaPz_in_out_vs_eta", 100, -3.0, 3.0, -100., 100.);
0143
0144   // Histogram to show the delta_Pl Out-In vs Z. The Z is the one
0145   // computed at the outermost hit.
0146   deltaPl_in_out_vs_z_ = ibook.bookProfile("DeltaPz_in_out_vs_z", "DeltaPz_in_out_vs_z", 600, -300, 300, -100., 100.);
0147
0148   // Histogram to show the delta_Pt Out-In in the Z-R plane. Z and R
0149   // are related to the outermost hit.
0150   deltaPt_in_out_2d_ = ibook.bookProfile2D("DeltaPt 2D", "DeltaPt 2D", 600, -300., 300, 120, 0., 120., -100., 100.);
0151
0152   // Histogram to show the distribution of p vs eta for all tracks.
0153   P_vs_eta_2d_ = ibook.book2D("P_vs_eta_2d", "P_vs_eta_2d", 100, -3.0, 3.0, 100., 0., 5.);
0154   char title[50];
0155   char key[20];
0156   for (unsigned int det = 1; det < sDETS.size(); ++det) {
0157     for (unsigned int sub_det = 1; sub_det <= trackerGeometry.numberOfLayers(det); ++sub_det) {
0158       memset(title, 0, sizeof(title));
0159       snprintf(title, sizeof(title), "Original_RadLen_vs_Eta_%s%d", sDETS[det].data(), sub_det);
0160       snprintf(key, sizeof(key), "%s%d", sDETS[det].data(), sub_det);
0161       histosOriEta_.insert(
0162           make_pair<string, MonitorElement *>(key, ibook.bookProfile(title, title, 250, -5.0, 5.0, 0., 1.)));
0163       snprintf(title, sizeof(title), "RadLen_vs_Eta_%s%d", sDETS[det].data(), sub_det);
0164       histosEta_.insert(
0165           make_pair<string, MonitorElement *>(key, ibook.bookProfile(title, title, 250, -5.0, 5.0, 0., 1.)));
0166     }
0167   }
0168 }
0169
0170 //-------------------------------------------------------------------------
0171 void TrackingRecoMaterialAnalyser::analyze(const edm::Event &event, const edm::EventSetup &setup) {
0172   using namespace edm;
0173   using namespace reco;
0174   using namespace std;
0175
0176   refitter_.setServices(setup);
0177
0178   Handle<TrackCollection> tracks = event.getHandle(tracksToken_);
0179   Handle<VertexCollection> vertices = event.getHandle(verticesToken_);
0180
0181   // Get the TrackerTopology
0182   const TrackerTopology *const tTopo = &setup.getData(tTopoToken_);
0183
0184   // Get Tracks
0185   if (!tracks.isValid() || tracks->empty()) {
0186     LogInfo("TrackingRecoMaterialAnalyser") << "Invalid or empty track collection" << endl;
0187     return;
0188   }
0189
0190   // Track Selector
0191   auto selector = [](const Track &track, const reco::Vertex::Point &pv) -> bool {
0192     return (track.quality(track.qualityByName("highPurity")) && track.dxy(pv) < 0.01 &&
0193             track.hitPattern().numberOfLostTrackerHits(HitPattern::MISSING_OUTER_HITS) == 0);
0194   };
0195
0196   // Get BeamSpot
0197   Handle<BeamSpot> beamSpot = event.getHandle(beamspotToken_);
0198   // Bail out if missing
0199   if (!beamSpot.isValid())
0200     return;
0201
0202   reco::Vertex::Point pv(beamSpot->position());
0203   if (usePV_) {
0204     if (vertices.isValid() && !vertices->empty()) {
0205       // Since we need to use eta and Z information from the tracks, in
0206       // order not to have the reco material distribution washed out due
0207       // to geometrical effects, we need to confine the PV to some small
0208       // region.
0209       const Vertex &v = (*vertices)[0];
0210       if (!v.isFake() && v.ndof() > 4 && std::fabs(v.z()) < 24 && std::fabs(v.position().rho()) < 2)
0211         pv = v.position();
0212     }
0213   }
0214
0215   // Main idea:
0216   // * select first good tracks in input, according to reasonable criteria
0217   // * refit the tracks so that we have access to the TrajectoryMeasurements
0218   //   that internally have all the information about the TSOS on all
0219   //   crossed layers. We need the refit track and not the original one so
0220   //   that we are able to correctly compute the path travelled by the track
0221   //   within the detector, using its updated TSOS. The material description
0222   //   can in principle be derived also directly from the rechits, via the
0223   //   det()[(->)GeomDet *]->mediumProperties chain, but that would simply give the
0224   //   face values, not the "real" ones used while propagating the track.
0225   // * Loop on all measurements, extract the information about the TSOS,
0226   //   its surface and its mediumProperties
0227   // * Make plots for the untouched material properties, but also for the
0228   //   ones corrected by the track direction, since the material properties,
0229   //   according to the documentation, should refer to normal incidence of
0230   //   the track, which is seldom the case, according to the current direction
0231   TrajectoryStateOnSurface current_tsos;
0232   DetId current_det;
0233   for (auto const &track : *tracks) {
0234     if (!selector(track, pv))
0235       continue;
0236     auto const &inner = track.innerMomentum();
0237     auto const &outer = track.outerMomentum();
0238     deltaP_in_out_vs_eta_->Fill(inner.eta(), inner.R() - outer.R());
0239     deltaP_in_out_vs_z_->Fill(track.outerZ(), inner.R() - outer.R());
0240     deltaP_in_out_vs_eta_2d_->Fill(inner.eta(), inner.R() - outer.R());
0241     deltaP_in_out_vs_eta_vs_phi_2d_->Fill(inner.eta(), inner.phi(), inner.R() - outer.R());
0242     deltaP_in_out_vs_z_2d_->Fill(track.outerZ(), inner.R() - outer.R());
0243     deltaPt_in_out_vs_eta_->Fill(inner.eta(), inner.rho() - outer.rho());
0244     deltaPt_in_out_vs_z_->Fill(track.outerZ(), inner.rho() - outer.rho());
0245     deltaPl_in_out_vs_eta_->Fill(inner.eta(), inner.z() - outer.z());
0246     deltaPl_in_out_vs_z_->Fill(track.outerZ(), inner.z() - outer.z());
0247     deltaPt_in_out_2d_->Fill(track.outerZ(), track.outerPosition().rho(), inner.rho() - outer.rho());
0248     P_vs_eta_2d_->Fill(track.eta(), track.p());
0249     vector<Trajectory> traj = refitter_.transform(track);
0250     if (traj.size() > 1 || traj.empty())
0251       continue;
0252     for (auto const &tm : traj.front().measurements()) {
0253       if (tm.recHit().get() &&
0254           (tm.recHitR().type() == TrackingRecHit::valid || tm.recHitR().type() == TrackingRecHit::missing)) {
0255         current_tsos = tm.updatedState().isValid() ? tm.updatedState() : tm.forwardPredictedState();
0256         auto const &localP = current_tsos.localMomentum();
0257         current_det = tm.recHit()->geographicalId();
0258         const Surface &surface = current_tsos.surface();
0259         assert(tm.recHit()->surface() == &surface);
0260         if (!surface.mediumProperties().isValid()) {
0261           LogError("TrackingRecoMaterialAnalyser")
0262               << "Medium properties for material linked to detector"
0263               << " are invalid at: " << current_tsos.globalPosition() << " " << (SiStripDetId)current_det << endl;
0264           assert(0);
0265           continue;
0266         }
0267         float p2 = localP.mag2();
0268         float xf = std::abs(std::sqrt(p2) / localP.z());
0269         float ori_xi = surface.mediumProperties().xi();
0270         float ori_radLen = surface.mediumProperties().radLen();
0271         float xi = ori_xi * xf;
0272         float radLen = ori_radLen * xf;
0273
0274         // NOTA BENE: THIS ASSUMES THAT THE TRACKS HAVE BEEN PRODUCED
0275         // WITH SPLITTING, AS IS THE CASE FOR THE generalTracks
0276         // collection.
0277
0278         // Since there are double-sided (glued) modules all over the
0279         // tracker, the material budget has been internally
0280         // partitioned in two equal components, so that each single
0281         // layer will receive half of the correct radLen. For this
0282         // reason, only for the double-sided components, we rescale
0283         // the obtained radLen by 2.
0284
0285         // In particular see code here:
0286         // http://cmslxr.fnal.gov/dxr/CMSSW_8_0_5/source/Geometry/TrackerGeometryBuilder/src/TrackerGeomBuilderFromGeometricDet.cc#213
0287         // where, in the SiStrip Tracker, if the module has a partner
0288         // (i.e. it's a glued detector) the plane is built with a
0289         // scaling of 0.5. The actual plane is built few lines below:
0290         // http://cmslxr.fnal.gov/dxr/CMSSW_8_0_5/source/Geometry/TrackerGeometryBuilder/src/TrackerGeomBuilderFromGeometricDet.cc#287
0291
0292         if (isDoubleSided(tTopo, current_det)) {
0293           LogTrace("TrackingRecoMaterialAnalyser")
0294               << "Eta: " << track.eta() << " " << sDETS[current_det.subdetId()] << tTopo->layer(current_det)
0295               << " has ori_radLen: " << ori_radLen << " and ori_xi: " << xi << " and has radLen: " << radLen
0296               << "  and xi: " << xi << endl;
0297           ori_radLen *= 2.;
0298           radLen *= 2.;
0299         }
0300
0301         histosOriEta_[sDETS[current_det.subdetId()] + to_string(tTopo->layer(current_det))]->Fill(
0302             current_tsos.globalPosition().eta(), ori_radLen);
0303         histosEta_[sDETS[current_det.subdetId()] + to_string(tTopo->layer(current_det))]->Fill(
0304             current_tsos.globalPosition().eta(), radLen);
0305         histo_RZ_Ori_->Fill(current_tsos.globalPosition().z(), current_tsos.globalPosition().perp(), ori_radLen);
0306         histo_RZ_->Fill(current_tsos.globalPosition().z(), current_tsos.globalPosition().perp(), radLen);
0307         LogInfo("TrackingRecoMaterialAnalyser")
0308             << "Eta: " << track.eta() << " " << sDETS[current_det.subdetId()] << tTopo->layer(current_det)
0309             << " has ori_radLen: " << ori_radLen << " and ori_xi: " << xi << " and has radLen: " << radLen
0310             << "  and xi: " << xi << endl;
0311       }
0312     }
0313   }
0314 }
0315
0316 //-------------------------------------------------------------------------
0317 // define as a plugin
0318 #include "FWCore/PluginManager/interface/ModuleDef.h"
0319 #include "FWCore/Framework/interface/MakerMacros.h"
0320 DEFINE_FWK_MODULE(TrackingRecoMaterialAnalyser);