File indexing completed on 2022-02-14 23:39:49
0001
0002
0003 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/ConsumesCollector.h"
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Utilities/interface/EDGetToken.h"
0010
0011 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0012 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0013
0014 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0015 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
0016 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
0017 #include "DataFormats/SiPixelCluster/interface/SiPixelClusterShapeCache.h"
0018
0019 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0020 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0021
0022 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0023 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0024
0025 #include "RecoPixelVertexing/PixelLowPtUtilities/interface/ClusterShapeHitFilter.h"
0026
0027 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0028
0029 #include "TROOT.h"
0030 #include "TFile.h"
0031 #include "TH1F.h"
0032 #include "TH2F.h"
0033
0034 #include <utility>
0035 #include <vector>
0036 #include <fstream>
0037 #include <memory>
0038
0039 using namespace std;
0040
0041
0042 #define exMax 10
0043 #define eyMax 15
0044
0045
0046 #define ewMax 40
0047
0048
0049 class ClusterShapeExtractor : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
0050 public:
0051 explicit ClusterShapeExtractor(const edm::ParameterSet& pset);
0052 ~ClusterShapeExtractor();
0053 virtual void beginRun(const edm::Run& run, const edm::EventSetup& es) override;
0054 virtual void analyze(const edm::Event& ev, const edm::EventSetup& es) override;
0055 virtual void endRun(const edm::Run& run, const edm::EventSetup& es) override{};
0056 virtual void endJob() override;
0057
0058 private:
0059 bool isSuitable(const PSimHit& simHit);
0060
0061
0062 void processSim(const SiPixelRecHit& recHit,
0063 const PSimHit& simHit,
0064 const SiPixelClusterShapeCache& clusterShapeCache,
0065 vector<TH2F*>& hspc);
0066 void processSim(const SiStripRecHit2D& recHit, const PSimHit& simHit, vector<TH1F*>& hssc);
0067
0068
0069 void processRec(const SiPixelRecHit& recHit,
0070 LocalVector ldir,
0071 const SiPixelClusterShapeCache& clusterShapeCache,
0072 vector<TH2F*>& hrpc);
0073 void processRec(const SiStripRecHit2D& recHit, LocalVector ldir, vector<TH1F*>& hrsc);
0074
0075 bool checkSimHits(const TrackingRecHit& recHit, PSimHit& simHit, pair<unsigned int, float>& key);
0076
0077 void processPixelRecHits(const SiPixelRecHitCollection::DataContainer* recHits,
0078 const SiPixelClusterShapeCache& clusterShapeCache);
0079 void processStripRecHits(const SiStripRecHit2DCollection::DataContainer* recHits);
0080 void processMatchedRecHits(const SiStripMatchedRecHit2DCollection::DataContainer* recHits);
0081
0082 void analyzeSimHits(const edm::Event& ev, const edm::EventSetup& es);
0083 void analyzeRecTracks(const edm::Event& ev, const edm::EventSetup& es);
0084
0085
0086 edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken;
0087 edm::ESGetToken<ClusterShapeHitFilter, CkfComponentsRecord> clusterShapeToken;
0088
0089 TFile* file;
0090
0091 string trackProducer;
0092 bool hasSimHits;
0093 bool hasRecTracks;
0094
0095 edm::EDGetTokenT<edmNew::DetSetVector<SiPixelRecHit>> pixelRecHits_token;
0096 edm::EDGetTokenT<SiPixelClusterShapeCache> clusterShapeCache_token;
0097
0098 const TrackerGeometry* theTracker;
0099 std::unique_ptr<TrackerHitAssociator> theHitAssociator;
0100 TrackerHitAssociator::Config trackerHitAssociatorConfig_;
0101 const ClusterShapeHitFilter* theClusterShape;
0102
0103 vector<TH2F*> hspc;
0104 vector<TH1F*> hssc;
0105
0106 vector<TH2F*> hrpc;
0107 vector<TH1F*> hrsc;
0108 };
0109
0110
0111 void ClusterShapeExtractor::beginRun(const edm::Run& run, const edm::EventSetup& es) {
0112
0113 theTracker = &es.getData(geomToken);
0114
0115
0116 theClusterShape = &es.getData(clusterShapeToken);
0117
0118
0119 char histName[256];
0120
0121
0122 for (int subdet = 0; subdet <= 1; subdet++) {
0123 for (int ex = 0; ex <= exMax; ex++)
0124 for (int ey = 0; ey <= eyMax; ey++) {
0125 sprintf(histName, "hspc_%d_%d_%d", subdet, ex, ey);
0126 hspc.push_back(new TH2F(histName,
0127 histName,
0128 10 * 2 * (exMax + 2),
0129 -(exMax + 2),
0130 (exMax + 2),
0131 10 * 2 * (eyMax + 2),
0132 -(eyMax + 2),
0133 (eyMax + 2)));
0134
0135 sprintf(histName, "hrpc_%d_%d_%d", subdet, ex, ey);
0136 hrpc.push_back(new TH2F(histName,
0137 histName,
0138 10 * 2 * (exMax + 2),
0139 -(exMax + 2),
0140 (exMax + 2),
0141 10 * 2 * (eyMax + 2),
0142 -(eyMax + 2),
0143 (eyMax + 2)));
0144 }
0145 }
0146
0147
0148 for (int ew = 0; ew <= ewMax; ew++) {
0149 sprintf(histName, "hssc_%d", ew);
0150 hssc.push_back(new TH1F(histName, histName, 10 * 2 * (ewMax * 2), -(ewMax * 2), (ewMax * 2)));
0151
0152 sprintf(histName, "hrsc_%d", ew);
0153 hrsc.push_back(new TH1F(histName, histName, 10 * 2 * (ewMax * 2), -(ewMax * 2), (ewMax * 2)));
0154 }
0155 }
0156
0157
0158 ClusterShapeExtractor::ClusterShapeExtractor(const edm::ParameterSet& pset)
0159 : geomToken(esConsumes<edm::Transition::BeginRun>()),
0160 clusterShapeToken(esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", "ClusterShapeHitFilter"))),
0161 pixelRecHits_token(consumes<edmNew::DetSetVector<SiPixelRecHit>>(edm::InputTag("siPixelRecHits"))),
0162 clusterShapeCache_token(
0163 consumes<SiPixelClusterShapeCache>(pset.getParameter<edm::InputTag>("clusterShapeCacheSrc"))),
0164 trackerHitAssociatorConfig_(pset, consumesCollector()) {
0165 trackProducer = pset.getParameter<string>("trackProducer");
0166 hasSimHits = pset.getParameter<bool>("hasSimHits");
0167 hasRecTracks = pset.getParameter<bool>("hasRecTracks");
0168
0169 file = new TFile("clusterShape.root", "RECREATE");
0170 file->cd();
0171 }
0172
0173
0174 void ClusterShapeExtractor::endJob() {
0175 typedef vector<TH2F*>::const_iterator H2I;
0176 typedef vector<TH1F*>::const_iterator H1I;
0177
0178 file->cd();
0179
0180
0181 for (H2I h = hspc.begin(); h != hspc.end(); h++)
0182 (*h)->Write();
0183 for (H1I h = hssc.begin(); h != hssc.end(); h++)
0184 (*h)->Write();
0185
0186
0187 for (H2I h = hrpc.begin(); h != hrpc.end(); h++)
0188 (*h)->Write();
0189 for (H1I h = hrsc.begin(); h != hrsc.end(); h++)
0190 (*h)->Write();
0191
0192 file->Close();
0193 }
0194
0195
0196 ClusterShapeExtractor::~ClusterShapeExtractor() {
0197
0198 }
0199
0200
0201 bool ClusterShapeExtractor::isSuitable(const PSimHit& simHit) {
0202
0203 DetId id = DetId(simHit.detUnitId());
0204 const GeomDetUnit* gdu = theTracker->idToDetUnit(id);
0205 if (gdu == 0)
0206 throw cms::Exception("MissingData") << "Missing DetUnit for detid " << id.rawId() << "\n" << std::endl;
0207 GlobalVector gvec = theTracker->idToDetUnit(id)->position() - GlobalPoint(0, 0, 0);
0208 LocalVector lvec = theTracker->idToDetUnit(id)->toLocal(gvec);
0209 LocalVector ldir = simHit.exitPoint() - simHit.entryPoint();
0210
0211 bool isOutgoing = (lvec.z() * ldir.z() > 0);
0212
0213 static const std::set<unsigned> RelevantProcesses = {0};
0214
0215 const bool isRelevant = RelevantProcesses.count(simHit.processType());
0216
0217
0218
0219
0220
0221 bool isFast = (simHit.momentumAtEntry().perp() > 0.050);
0222
0223
0224 return (isOutgoing && isRelevant && isFast);
0225 }
0226
0227
0228 void ClusterShapeExtractor::processRec(const SiStripRecHit2D& recHit, LocalVector ldir, vector<TH1F*>& histo) {
0229 int meas;
0230 float pred;
0231
0232 if (theClusterShape->getSizes(recHit, LocalPoint(0, 0, 0), ldir, meas, pred))
0233 if (meas <= ewMax)
0234 histo[meas]->Fill(pred);
0235 }
0236
0237
0238 void ClusterShapeExtractor::processRec(const SiPixelRecHit& recHit,
0239 LocalVector ldir,
0240 const SiPixelClusterShapeCache& clusterShapeCache,
0241 vector<TH2F*>& histo) {
0242 int part;
0243 ClusterData::ArrayType meas;
0244 pair<float, float> pred;
0245
0246 if (theClusterShape->getSizes(recHit, ldir, clusterShapeCache, part, meas, pred))
0247 if (meas.size() == 1)
0248 if (meas.front().first <= exMax && meas.front().second <= eyMax) {
0249 int i = (part * (exMax + 1) + meas.front().first) * (eyMax + 1) + meas.front().second;
0250 histo[i]->Fill(pred.first, pred.second);
0251 }
0252 }
0253
0254
0255 void ClusterShapeExtractor::processSim(const SiPixelRecHit& recHit,
0256 const PSimHit& simHit,
0257 const SiPixelClusterShapeCache& clusterShapeCache,
0258 vector<TH2F*>& histo) {
0259 LocalVector ldir = simHit.exitPoint() - simHit.entryPoint();
0260 processRec(recHit, ldir, clusterShapeCache, histo);
0261 }
0262
0263
0264 void ClusterShapeExtractor::processSim(const SiStripRecHit2D& recHit, const PSimHit& simHit, vector<TH1F*>& histo) {
0265 LocalVector ldir = simHit.exitPoint() - simHit.entryPoint();
0266 processRec(recHit, ldir, histo);
0267 }
0268
0269
0270 bool ClusterShapeExtractor::checkSimHits(const TrackingRecHit& recHit,
0271 PSimHit& simHit,
0272 pair<unsigned int, float>& key) {
0273 vector<PSimHit> simHits = theHitAssociator->associateHit(recHit);
0274
0275
0276 if (simHits.size() == 1) {
0277 simHit = simHits[0];
0278
0279 if (isSuitable(simHit)) {
0280 key = pair<unsigned int, float>(simHit.trackId(), simHit.timeOfFlight());
0281 return true;
0282 }
0283 }
0284
0285 return false;
0286 }
0287
0288
0289 void ClusterShapeExtractor::processPixelRecHits(const SiPixelRecHitCollection::DataContainer* recHits,
0290 const SiPixelClusterShapeCache& clusterShapeCache) {
0291 map<pair<unsigned int, float>, const SiPixelRecHit*> simHitMap;
0292
0293 PSimHit simHit;
0294 pair<unsigned int, float> key;
0295 size_t counter = 0, counter_2 = 0;
0296
0297 for (SiPixelRecHitCollection::DataContainer::const_iterator recHit = recHits->begin(); recHit != recHits->end();
0298 recHit++)
0299 if (checkSimHits(*recHit, simHit, key)) {
0300
0301 if (simHitMap.count(key) == 0) {
0302 simHitMap[key] = &(*recHit);
0303 } else if (recHit->cluster()->size() > simHitMap[key]->cluster()->size())
0304 simHitMap[key] = &(*recHit);
0305 ++counter_2;
0306 }
0307
0308 for (SiPixelRecHitCollection::DataContainer::const_iterator recHit = recHits->begin(); recHit != recHits->end();
0309 recHit++)
0310 if (checkSimHits(*recHit, simHit, key)) {
0311
0312 if (&(*recHit) == simHitMap[key]) {
0313 processSim(*recHit, simHit, clusterShapeCache, hspc);
0314 ++counter;
0315 }
0316 }
0317
0318
0319 }
0320
0321
0322 void ClusterShapeExtractor::processStripRecHits(const SiStripRecHit2DCollection::DataContainer* recHits) {
0323 map<pair<unsigned int, float>, const SiStripRecHit2D*> simHitMap;
0324
0325 PSimHit simHit;
0326 pair<unsigned int, float> key;
0327
0328 for (SiStripRecHit2DCollection::DataContainer::const_iterator recHit = recHits->begin(); recHit != recHits->end();
0329 recHit++)
0330 if (checkSimHits(*recHit, simHit, key)) {
0331
0332 if (simHitMap.count(key) == 0)
0333 simHitMap[key] = &(*recHit);
0334 else if (recHit->cluster()->amplitudes().size() > simHitMap[key]->cluster()->amplitudes().size())
0335 simHitMap[key] = &(*recHit);
0336 }
0337
0338 for (SiStripRecHit2DCollection::DataContainer::const_iterator recHit = recHits->begin(); recHit != recHits->end();
0339 recHit++)
0340 if (checkSimHits(*recHit, simHit, key)) {
0341
0342 if (&(*recHit) == simHitMap[key])
0343 processSim(*recHit, simHit, hssc);
0344 }
0345 }
0346
0347
0348 void ClusterShapeExtractor::processMatchedRecHits(const SiStripMatchedRecHit2DCollection::DataContainer* recHits) {
0349 map<pair<unsigned int, float>, const SiStripRecHit2D*> simHitMap;
0350
0351 std::vector<SiStripRecHit2D> cache;
0352 cache.reserve(2 * recHits->size());
0353 PSimHit simHit;
0354 pair<unsigned int, float> key;
0355
0356 const SiStripRecHit2D* recHit;
0357
0358 for (SiStripMatchedRecHit2DCollection::DataContainer::const_iterator matched = recHits->begin();
0359 matched != recHits->end();
0360 matched++) {
0361 cache.push_back(matched->monoHit());
0362 recHit = &cache.back();
0363 if (checkSimHits(*recHit, simHit, key)) {
0364
0365 if (simHitMap.count(key) == 0)
0366 simHitMap[key] = &(*recHit);
0367 else if (recHit->cluster()->amplitudes().size() > simHitMap[key]->cluster()->amplitudes().size())
0368 simHitMap[key] = &(*recHit);
0369 }
0370 cache.push_back(matched->stereoHit());
0371 recHit = &cache.back();
0372 if (checkSimHits(*recHit, simHit, key)) {
0373
0374 if (simHitMap.count(key) == 0)
0375 simHitMap[key] = &(*recHit);
0376 else if (recHit->cluster()->amplitudes().size() > simHitMap[key]->cluster()->amplitudes().size())
0377 simHitMap[key] = &(*recHit);
0378 }
0379 }
0380
0381 for (SiStripMatchedRecHit2DCollection::DataContainer::const_iterator matched = recHits->begin();
0382 matched != recHits->end();
0383 matched++) {
0384 auto recHit = matched->monoHit();
0385 if (checkSimHits(recHit, simHit, key)) {
0386
0387 if (recHit.omniCluster() == simHitMap[key]->omniCluster())
0388 processSim(recHit, simHit, hssc);
0389 }
0390
0391 recHit = matched->stereoHit();
0392 if (checkSimHits(recHit, simHit, key)) {
0393
0394 if (recHit.omniCluster() == simHitMap[key]->omniCluster())
0395 processSim(recHit, simHit, hssc);
0396 }
0397 }
0398 }
0399
0400
0401 void ClusterShapeExtractor::analyzeSimHits(const edm::Event& ev, const edm::EventSetup& es) {
0402
0403 theHitAssociator.reset(new TrackerHitAssociator(ev, trackerHitAssociatorConfig_));
0404
0405
0406 {
0407 edm::Handle<SiPixelRecHitCollection> coll;
0408
0409 ev.getByToken(pixelRecHits_token, coll);
0410
0411 edm::Handle<SiPixelClusterShapeCache> clusterShapeCache;
0412 ev.getByToken(clusterShapeCache_token, clusterShapeCache);
0413
0414 const SiPixelRecHitCollection::DataContainer* recHits = &coll.product()->data();
0415 processPixelRecHits(recHits, *clusterShapeCache);
0416 }
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439
0440
0441
0442
0443
0444
0445
0446
0447 }
0448
0449
0450 void ClusterShapeExtractor::analyzeRecTracks(const edm::Event& ev, const edm::EventSetup& es) {
0451
0452 edm::Handle<vector<Trajectory>> trajeHandle;
0453 ev.getByLabel(trackProducer, trajeHandle);
0454 const vector<Trajectory>& trajeCollection = *(trajeHandle.product());
0455
0456 edm::Handle<SiPixelClusterShapeCache> clusterShapeCache;
0457 ev.getByToken(clusterShapeCache_token, clusterShapeCache);
0458
0459
0460 for (vector<Trajectory>::const_iterator trajectory = trajeCollection.begin(); trajectory != trajeCollection.end();
0461 trajectory++)
0462 for (vector<TrajectoryMeasurement>::const_iterator meas = trajectory->measurements().begin();
0463 meas != trajectory->measurements().end();
0464 meas++) {
0465 const TrackingRecHit* recHit = meas->recHit()->hit();
0466 DetId id = recHit->geographicalId();
0467
0468 if (recHit->isValid()) {
0469 LocalVector ldir = meas->updatedState().localDirection();
0470
0471 if (GeomDetEnumerators::isTrackerPixel(theTracker->geomDetSubDetector(id.subdetId()))) {
0472
0473 const SiPixelRecHit* pixelRecHit = dynamic_cast<const SiPixelRecHit*>(recHit);
0474
0475 if (pixelRecHit != 0)
0476 processRec(*pixelRecHit, ldir, *clusterShapeCache, hrpc);
0477 } else if (GeomDetEnumerators::isTrackerStrip(theTracker->geomDetSubDetector(id.subdetId()))) {
0478
0479 const SiStripMatchedRecHit2D* stripMatchedRecHit = dynamic_cast<const SiStripMatchedRecHit2D*>(recHit);
0480 const ProjectedSiStripRecHit2D* stripProjectedRecHit = dynamic_cast<const ProjectedSiStripRecHit2D*>(recHit);
0481 const SiStripRecHit2D* stripRecHit = dynamic_cast<const SiStripRecHit2D*>(recHit);
0482
0483 if (stripMatchedRecHit != 0) {
0484 processRec(stripMatchedRecHit->monoHit(), ldir, hrsc);
0485 processRec(stripMatchedRecHit->stereoHit(), ldir, hrsc);
0486 }
0487
0488 if (stripProjectedRecHit != 0)
0489 processRec(stripProjectedRecHit->originalHit(), ldir, hrsc);
0490
0491 if (stripRecHit != 0)
0492 processRec(*stripRecHit, ldir, hrsc);
0493 }
0494 }
0495 }
0496 }
0497
0498
0499 void ClusterShapeExtractor::analyze(const edm::Event& ev, const edm::EventSetup& es) {
0500 if (hasSimHits) {
0501 LogTrace("MinBiasTracking") << " [ClusterShape] analyze simHits, recHits";
0502 analyzeSimHits(ev, es);
0503 }
0504
0505 if (hasRecTracks) {
0506 LogTrace("MinBiasTracking") << " [ClusterShape] analyze recHits on recTracks";
0507 analyzeRecTracks(ev, es);
0508 }
0509 }
0510
0511 DEFINE_FWK_MODULE(ClusterShapeExtractor);