File indexing completed on 2024-04-06 12:09:09
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
0029
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
0103
0104
0105 histo_RZ_Ori_ = ibook.bookProfile2D("OriRadLen", "Original_RadLen", 600, -300., 300, 120, 0., 120., 0., 1.);
0106
0107
0108
0109
0110 histo_RZ_ = ibook.bookProfile2D("RadLen", "RadLen", 600, -300., 300, 120, 0., 120., 0., 1.);
0111
0112
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
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
0121
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
0125
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
0131
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
0136
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
0140
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
0145
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
0149
0150 deltaPt_in_out_2d_ = ibook.bookProfile2D("DeltaPt 2D", "DeltaPt 2D", 600, -300., 300, 120, 0., 120., -100., 100.);
0151
0152
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
0182 const TrackerTopology *const tTopo = &setup.getData(tTopoToken_);
0183
0184
0185 if (!tracks.isValid() || tracks->empty()) {
0186 LogInfo("TrackingRecoMaterialAnalyser") << "Invalid or empty track collection" << endl;
0187 return;
0188 }
0189
0190
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
0197 Handle<BeamSpot> beamSpot = event.getHandle(beamspotToken_);
0198
0199 if (!beamSpot.isValid())
0200 return;
0201
0202 reco::Vertex::Point pv(beamSpot->position());
0203 if (usePV_) {
0204 if (vertices.isValid() && !vertices->empty()) {
0205
0206
0207
0208
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
0216
0217
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230
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
0275
0276
0277
0278
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288
0289
0290
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
0318 #include "FWCore/PluginManager/interface/ModuleDef.h"
0319 #include "FWCore/Framework/interface/MakerMacros.h"
0320 DEFINE_FWK_MODULE(TrackingRecoMaterialAnalyser);