File indexing completed on 2024-09-07 04:38:00
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 "RecoTracker/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
0296 for (SiPixelRecHitCollection::DataContainer::const_iterator recHit = recHits->begin(); recHit != recHits->end();
0297 recHit++)
0298 if (checkSimHits(*recHit, simHit, key)) {
0299
0300 if (simHitMap.count(key) == 0) {
0301 simHitMap[key] = &(*recHit);
0302 } else if (recHit->cluster()->size() > simHitMap[key]->cluster()->size())
0303 simHitMap[key] = &(*recHit);
0304 }
0305
0306 for (SiPixelRecHitCollection::DataContainer::const_iterator recHit = recHits->begin(); recHit != recHits->end();
0307 recHit++)
0308 if (checkSimHits(*recHit, simHit, key)) {
0309
0310 if (&(*recHit) == simHitMap[key]) {
0311 processSim(*recHit, simHit, clusterShapeCache, hspc);
0312 }
0313 }
0314 }
0315
0316
0317 void ClusterShapeExtractor::processStripRecHits(const SiStripRecHit2DCollection::DataContainer* recHits) {
0318 map<pair<unsigned int, float>, const SiStripRecHit2D*> simHitMap;
0319
0320 PSimHit simHit;
0321 pair<unsigned int, float> key;
0322
0323 for (SiStripRecHit2DCollection::DataContainer::const_iterator recHit = recHits->begin(); recHit != recHits->end();
0324 recHit++)
0325 if (checkSimHits(*recHit, simHit, key)) {
0326
0327 if (simHitMap.count(key) == 0)
0328 simHitMap[key] = &(*recHit);
0329 else if (recHit->cluster()->amplitudes().size() > simHitMap[key]->cluster()->amplitudes().size())
0330 simHitMap[key] = &(*recHit);
0331 }
0332
0333 for (SiStripRecHit2DCollection::DataContainer::const_iterator recHit = recHits->begin(); recHit != recHits->end();
0334 recHit++)
0335 if (checkSimHits(*recHit, simHit, key)) {
0336
0337 if (&(*recHit) == simHitMap[key])
0338 processSim(*recHit, simHit, hssc);
0339 }
0340 }
0341
0342
0343 void ClusterShapeExtractor::processMatchedRecHits(const SiStripMatchedRecHit2DCollection::DataContainer* recHits) {
0344 map<pair<unsigned int, float>, const SiStripRecHit2D*> simHitMap;
0345
0346 std::vector<SiStripRecHit2D> cache;
0347 cache.reserve(2 * recHits->size());
0348 PSimHit simHit;
0349 pair<unsigned int, float> key;
0350
0351 const SiStripRecHit2D* recHit;
0352
0353 for (SiStripMatchedRecHit2DCollection::DataContainer::const_iterator matched = recHits->begin();
0354 matched != recHits->end();
0355 matched++) {
0356 cache.push_back(matched->monoHit());
0357 recHit = &cache.back();
0358 if (checkSimHits(*recHit, simHit, key)) {
0359
0360 if (simHitMap.count(key) == 0)
0361 simHitMap[key] = &(*recHit);
0362 else if (recHit->cluster()->amplitudes().size() > simHitMap[key]->cluster()->amplitudes().size())
0363 simHitMap[key] = &(*recHit);
0364 }
0365 cache.push_back(matched->stereoHit());
0366 recHit = &cache.back();
0367 if (checkSimHits(*recHit, simHit, key)) {
0368
0369 if (simHitMap.count(key) == 0)
0370 simHitMap[key] = &(*recHit);
0371 else if (recHit->cluster()->amplitudes().size() > simHitMap[key]->cluster()->amplitudes().size())
0372 simHitMap[key] = &(*recHit);
0373 }
0374 }
0375
0376 for (SiStripMatchedRecHit2DCollection::DataContainer::const_iterator matched = recHits->begin();
0377 matched != recHits->end();
0378 matched++) {
0379 auto recHit = matched->monoHit();
0380 if (checkSimHits(recHit, simHit, key)) {
0381
0382 if (recHit.omniCluster() == simHitMap[key]->omniCluster())
0383 processSim(recHit, simHit, hssc);
0384 }
0385
0386 recHit = matched->stereoHit();
0387 if (checkSimHits(recHit, simHit, key)) {
0388
0389 if (recHit.omniCluster() == simHitMap[key]->omniCluster())
0390 processSim(recHit, simHit, hssc);
0391 }
0392 }
0393 }
0394
0395
0396 void ClusterShapeExtractor::analyzeSimHits(const edm::Event& ev, const edm::EventSetup& es) {
0397
0398 theHitAssociator.reset(new TrackerHitAssociator(ev, trackerHitAssociatorConfig_));
0399
0400
0401 {
0402 edm::Handle<SiPixelRecHitCollection> coll;
0403
0404 ev.getByToken(pixelRecHits_token, coll);
0405
0406 edm::Handle<SiPixelClusterShapeCache> clusterShapeCache;
0407 ev.getByToken(clusterShapeCache_token, clusterShapeCache);
0408
0409 const SiPixelRecHitCollection::DataContainer* recHits = &coll.product()->data();
0410 processPixelRecHits(recHits, *clusterShapeCache);
0411 }
0412
0413
0414
0415
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 void ClusterShapeExtractor::analyzeRecTracks(const edm::Event& ev, const edm::EventSetup& es) {
0446
0447 edm::Handle<vector<Trajectory>> trajeHandle;
0448 ev.getByLabel(trackProducer, trajeHandle);
0449 const vector<Trajectory>& trajeCollection = *(trajeHandle.product());
0450
0451 edm::Handle<SiPixelClusterShapeCache> clusterShapeCache;
0452 ev.getByToken(clusterShapeCache_token, clusterShapeCache);
0453
0454
0455 for (vector<Trajectory>::const_iterator trajectory = trajeCollection.begin(); trajectory != trajeCollection.end();
0456 trajectory++)
0457 for (vector<TrajectoryMeasurement>::const_iterator meas = trajectory->measurements().begin();
0458 meas != trajectory->measurements().end();
0459 meas++) {
0460 const TrackingRecHit* recHit = meas->recHit()->hit();
0461 DetId id = recHit->geographicalId();
0462
0463 if (recHit->isValid()) {
0464 LocalVector ldir = meas->updatedState().localDirection();
0465
0466 if (GeomDetEnumerators::isTrackerPixel(theTracker->geomDetSubDetector(id.subdetId()))) {
0467
0468 const SiPixelRecHit* pixelRecHit = dynamic_cast<const SiPixelRecHit*>(recHit);
0469
0470 if (pixelRecHit != 0)
0471 processRec(*pixelRecHit, ldir, *clusterShapeCache, hrpc);
0472 } else if (GeomDetEnumerators::isTrackerStrip(theTracker->geomDetSubDetector(id.subdetId()))) {
0473
0474 const SiStripMatchedRecHit2D* stripMatchedRecHit = dynamic_cast<const SiStripMatchedRecHit2D*>(recHit);
0475 const ProjectedSiStripRecHit2D* stripProjectedRecHit = dynamic_cast<const ProjectedSiStripRecHit2D*>(recHit);
0476 const SiStripRecHit2D* stripRecHit = dynamic_cast<const SiStripRecHit2D*>(recHit);
0477
0478 if (stripMatchedRecHit != 0) {
0479 processRec(stripMatchedRecHit->monoHit(), ldir, hrsc);
0480 processRec(stripMatchedRecHit->stereoHit(), ldir, hrsc);
0481 }
0482
0483 if (stripProjectedRecHit != 0)
0484 processRec(stripProjectedRecHit->originalHit(), ldir, hrsc);
0485
0486 if (stripRecHit != 0)
0487 processRec(*stripRecHit, ldir, hrsc);
0488 }
0489 }
0490 }
0491 }
0492
0493
0494 void ClusterShapeExtractor::analyze(const edm::Event& ev, const edm::EventSetup& es) {
0495 if (hasSimHits) {
0496 LogTrace("MinBiasTracking") << " [ClusterShape] analyze simHits, recHits";
0497 analyzeSimHits(ev, es);
0498 }
0499
0500 if (hasRecTracks) {
0501 LogTrace("MinBiasTracking") << " [ClusterShape] analyze recHits on recTracks";
0502 analyzeRecTracks(ev, es);
0503 }
0504 }
0505
0506 DEFINE_FWK_MODULE(ClusterShapeExtractor);