File indexing completed on 2023-03-31 23:02:08
0001
0002
0003 #include "FWCore/Framework/interface/global/EDAnalyzer.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/ESHandle.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/SiPixelCluster/interface/SiPixelClusterShapeCache.h"
0016
0017 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0018 #include "DataFormats/TrackReco/interface/Track.h"
0019
0020 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0021 #include "RecoTracker/PixelLowPtUtilities/interface/ClusterShapeHitFilter.h"
0022
0023 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0024 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0025 #
0026
0027 #include "TROOT.h"
0028 #include "TFile.h"
0029 #include "TH1F.h"
0030 #include "TH2F.h"
0031
0032 #include <map>
0033 #include <memory>
0034
0035 #include <fstream>
0036 #include <vector>
0037
0038 #include <mutex>
0039
0040 using namespace std;
0041
0042
0043 #define exMax 10
0044 #define eyMax 15
0045
0046 namespace {
0047
0048 const std::set<unsigned> RelevantProcesses = {0};
0049
0050
0051 }
0052
0053
0054 class PixelClusterShapeExtractor final : public edm::global::EDAnalyzer<> {
0055 public:
0056 explicit PixelClusterShapeExtractor(const edm::ParameterSet& pset);
0057 void analyze(edm::StreamID, const edm::Event& evt, const edm::EventSetup&) const override;
0058 void endJob() override;
0059
0060 private:
0061 void init();
0062
0063 bool isSuitable(const PSimHit& simHit, const GeomDetUnit& gdu) const;
0064
0065
0066 void processSim(const SiPixelRecHit& recHit,
0067 ClusterShapeHitFilter const& theClusterFilter,
0068 const PSimHit& simHit,
0069 const SiPixelClusterShapeCache& clusterShapeCache,
0070 const vector<TH2F*>& histo) const;
0071
0072
0073 void processRec(const SiPixelRecHit& recHit,
0074 ClusterShapeHitFilter const& theFilter,
0075 LocalVector ldir,
0076 const SiPixelClusterShapeCache& clusterShapeCache,
0077 const vector<TH2F*>& histo) const;
0078
0079 bool checkSimHits(const TrackingRecHit& recHit,
0080 TrackerHitAssociator const& theAssociator,
0081 PSimHit& simHit,
0082 pair<unsigned int, float>& key,
0083 unsigned int& ss) const;
0084
0085 void processPixelRecHits(SiPixelRecHitCollection::DataContainer const& recHits,
0086 TrackerHitAssociator const& theAssociator,
0087 ClusterShapeHitFilter const& theFilter,
0088 SiPixelClusterShapeCache const& clusterShapeCache,
0089 const TrackerTopology& tkTpl) const;
0090
0091 void analyzeSimHits(const edm::Event& ev, const edm::EventSetup& es) const;
0092 void analyzeRecTracks(const edm::Event& ev, const edm::EventSetup& es) const;
0093
0094
0095 const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0096 const edm::ESGetToken<ClusterShapeHitFilter, CkfComponentsRecord> csfToken_;
0097
0098 TFile* file;
0099
0100 const bool hasSimHits;
0101 const bool hasRecTracks;
0102 const bool noBPIX1;
0103
0104 const edm::EDGetTokenT<reco::TrackCollection> tracks_token;
0105 const edm::EDGetTokenT<edmNew::DetSetVector<SiPixelRecHit>> pixelRecHits_token;
0106 const edm::EDGetTokenT<SiPixelClusterShapeCache> clusterShapeCache_token;
0107 const TrackerHitAssociator::Config trackerHitAssociatorConfig_;
0108
0109 using Lock = std::unique_lock<std::mutex>;
0110 std::unique_ptr<std::mutex[]> theMutex;
0111 std::vector<TH2F*> hspc;
0112 std::vector<TH2F*> hrpc;
0113 };
0114
0115
0116 void PixelClusterShapeExtractor::init() {
0117
0118 char histName[256];
0119
0120
0121 for (int subdet = 0; subdet <= 1; subdet++) {
0122 for (int ex = 0; ex <= exMax; ex++)
0123 for (int ey = 0; ey <= eyMax; ey++) {
0124 sprintf(histName, "hspc_%d_%d_%d", subdet, ex, ey);
0125 hspc.push_back(new TH2F(histName,
0126 histName,
0127 10 * 2 * (exMax + 2),
0128 -(exMax + 2),
0129 (exMax + 2),
0130 10 * 2 * (eyMax + 2),
0131 -(eyMax + 2),
0132 (eyMax + 2)));
0133
0134 sprintf(histName, "hrpc_%d_%d_%d", subdet, ex, ey);
0135 hrpc.push_back(new TH2F(histName,
0136 histName,
0137 10 * 2 * (exMax + 2),
0138 -(exMax + 2),
0139 (exMax + 2),
0140 10 * 2 * (eyMax + 2),
0141 -(eyMax + 2),
0142 (eyMax + 2)));
0143 }
0144 }
0145 theMutex = std::make_unique<std::mutex[]>(hspc.size());
0146 }
0147
0148
0149 PixelClusterShapeExtractor::PixelClusterShapeExtractor(const edm::ParameterSet& pset)
0150 : topoToken_(esConsumes()),
0151 csfToken_(esConsumes(edm::ESInputTag("", "ClusterShapeHitFilter"))),
0152 hasSimHits(pset.getParameter<bool>("hasSimHits")),
0153 hasRecTracks(pset.getParameter<bool>("hasRecTracks")),
0154 noBPIX1(pset.getParameter<bool>("noBPIX1")),
0155 tracks_token(hasRecTracks ? consumes<reco::TrackCollection>(pset.getParameter<edm::InputTag>("tracks"))
0156 : edm::EDGetTokenT<reco::TrackCollection>()),
0157 pixelRecHits_token(consumes<edmNew::DetSetVector<SiPixelRecHit>>(edm::InputTag("siPixelRecHits"))),
0158 clusterShapeCache_token(
0159 consumes<SiPixelClusterShapeCache>(pset.getParameter<edm::InputTag>("clusterShapeCacheSrc"))),
0160 trackerHitAssociatorConfig_(pset, consumesCollector()) {
0161 file = new TFile("clusterShape.root", "RECREATE");
0162 file->cd();
0163 init();
0164 }
0165
0166
0167 void PixelClusterShapeExtractor::endJob() {
0168 file->cd();
0169
0170
0171 for (auto h = hspc.begin(); h != hspc.end(); h++)
0172 (*h)->Write();
0173
0174
0175 for (auto h = hrpc.begin(); h != hrpc.end(); h++)
0176 (*h)->Write();
0177
0178 file->Close();
0179 }
0180
0181
0182 bool PixelClusterShapeExtractor::isSuitable(const PSimHit& simHit, const GeomDetUnit& gdu) const {
0183
0184
0185 GlobalVector gvec = gdu.position() - GlobalPoint(0, 0, 0);
0186 LocalVector lvec = gdu.toLocal(gvec);
0187 LocalVector ldir = simHit.exitPoint() - simHit.entryPoint();
0188
0189
0190 bool isOutgoing = std::abs(ldir.z()) > 0.01f && (lvec.z() * ldir.z() > 0);
0191
0192
0193 const bool isRelevant = RelevantProcesses.count(simHit.processType());
0194
0195
0196
0197
0198 constexpr float ptCut2 = 0.2 * 0.2;
0199
0200 bool isFast = (simHit.momentumAtEntry().perp2() > ptCut2);
0201
0202
0203 return (isOutgoing && isRelevant && isFast);
0204 }
0205
0206
0207 void PixelClusterShapeExtractor::processRec(const SiPixelRecHit& recHit,
0208 ClusterShapeHitFilter const& theClusterShape,
0209 LocalVector ldir,
0210 const SiPixelClusterShapeCache& clusterShapeCache,
0211 const vector<TH2F*>& histo) const {
0212 int part;
0213 ClusterData::ArrayType meas;
0214 pair<float, float> pred;
0215
0216 if (theClusterShape.getSizes(recHit, ldir, clusterShapeCache, part, meas, pred))
0217 if (meas.size() == 1)
0218 if (meas.front().first <= exMax && meas.front().second <= eyMax) {
0219 int i = (part * (exMax + 1) + meas.front().first) * (eyMax + 1) + meas.front().second;
0220 #ifdef DO_DEBUG
0221 if (meas.front().second == 0 && std::abs(pred.second) > 3) {
0222 Lock lock(theMutex[0]);
0223 int id = recHit.geographicalId();
0224 std::cout << id << " bigpred " << meas.front().first << '/' << meas.front().second << ' ' << pred.first << '/'
0225 << pred.second << ' ' << ldir << ' ' << ldir.mag() << std::endl;
0226 }
0227 #endif
0228 Lock lock(theMutex[i]);
0229 histo[i]->Fill(pred.first, pred.second);
0230 }
0231 }
0232
0233
0234 void PixelClusterShapeExtractor::processSim(const SiPixelRecHit& recHit,
0235 ClusterShapeHitFilter const& theClusterFilter,
0236 const PSimHit& simHit,
0237 const SiPixelClusterShapeCache& clusterShapeCache,
0238 const vector<TH2F*>& histo) const {
0239 LocalVector ldir = simHit.exitPoint() - simHit.entryPoint();
0240 processRec(recHit, theClusterFilter, ldir, clusterShapeCache, histo);
0241 }
0242
0243
0244 bool PixelClusterShapeExtractor::checkSimHits(const TrackingRecHit& recHit,
0245 TrackerHitAssociator const& theHitAssociator,
0246 PSimHit& simHit,
0247 pair<unsigned int, float>& key,
0248 unsigned int& ss) const {
0249 auto const& simHits = theHitAssociator.associateHit(recHit);
0250
0251
0252 for (auto const& sh : simHits) {
0253 if (isSuitable(sh, *recHit.detUnit())) {
0254 simHit = sh;
0255 key = {simHit.trackId(), simHit.timeOfFlight()};
0256 ss = simHits.size();
0257 return true;
0258 }
0259 }
0260
0261 return false;
0262 }
0263
0264
0265 void PixelClusterShapeExtractor::processPixelRecHits(const SiPixelRecHitCollection::DataContainer& recHits,
0266 TrackerHitAssociator const& theHitAssociator,
0267 ClusterShapeHitFilter const& theFilter,
0268 const SiPixelClusterShapeCache& clusterShapeCache,
0269 const TrackerTopology& tkTpl) const {
0270 struct Elem {
0271 const SiPixelRecHit* rhit;
0272 PSimHit shit;
0273 unsigned int size;
0274 };
0275 std::map<pair<unsigned int, float>, Elem> simHitMap;
0276
0277 PSimHit simHit;
0278 pair<unsigned int, float> key;
0279 unsigned int ss;
0280
0281 for (auto const& recHit : recHits) {
0282 if (noBPIX1 && tkTpl.pxbLayer(recHit.geographicalId()) == 1)
0283 continue;
0284 if (!checkSimHits(recHit, theHitAssociator, simHit, key, ss))
0285 continue;
0286
0287 if (simHitMap.count(key) == 0) {
0288 simHitMap[key] = {&recHit, simHit, ss};
0289 } else if (recHit.cluster()->size() > simHitMap[key].rhit->cluster()->size())
0290 simHitMap[key] = {&recHit, simHit, std::max(ss, simHitMap[key].size)};
0291 }
0292 for (auto const& elem : simHitMap) {
0293
0294
0295
0296
0297
0298
0299
0300
0301
0302 if (elem.second.size == 1)
0303 processSim(*elem.second.rhit, theFilter, elem.second.shit, clusterShapeCache, hspc);
0304 }
0305 }
0306
0307
0308 void PixelClusterShapeExtractor::analyzeSimHits(const edm::Event& ev, const edm::EventSetup& es) const {
0309 auto const& theClusterShape = es.getData(csfToken_);
0310 auto const& tkTpl = es.getData(topoToken_);
0311
0312 edm::Handle<SiPixelClusterShapeCache> clusterShapeCache;
0313 ev.getByToken(clusterShapeCache_token, clusterShapeCache);
0314
0315
0316 auto theHitAssociator = std::make_unique<TrackerHitAssociator>(ev, trackerHitAssociatorConfig_);
0317
0318
0319 {
0320 edm::Handle<SiPixelRecHitCollection> coll;
0321 ev.getByToken(pixelRecHits_token, coll);
0322
0323 edm::Handle<SiPixelClusterShapeCache> clusterShapeCache;
0324 ev.getByToken(clusterShapeCache_token, clusterShapeCache);
0325
0326 auto const& recHits = coll.product()->data();
0327 processPixelRecHits(recHits, *theHitAssociator, theClusterShape, *clusterShapeCache, tkTpl);
0328 }
0329 }
0330
0331
0332 void PixelClusterShapeExtractor::analyzeRecTracks(const edm::Event& ev, const edm::EventSetup& es) const {
0333 auto const& theClusterShape = es.getData(csfToken_);
0334 auto const& tkTpl = es.getData(topoToken_);
0335
0336
0337 edm::Handle<reco::TrackCollection> tracks;
0338 ev.getByToken(tracks_token, tracks);
0339
0340 edm::Handle<SiPixelClusterShapeCache> clusterShapeCache;
0341 ev.getByToken(clusterShapeCache_token, clusterShapeCache);
0342
0343 for (auto const& track : *tracks) {
0344 if (!track.quality(reco::Track::highPurity))
0345 continue;
0346 if (track.numberOfValidHits() < 8)
0347 continue;
0348 auto const& trajParams = track.extra()->trajParams();
0349 assert(trajParams.size() == track.recHitsSize());
0350 auto hb = track.recHitsBegin();
0351 for (unsigned int h = 0; h < track.recHitsSize(); h++) {
0352 auto recHit = *(hb + h);
0353 if (!recHit->isValid())
0354 continue;
0355 auto id = recHit->geographicalId();
0356 if (noBPIX1 && tkTpl.pxbLayer(id) == 1)
0357 continue;
0358
0359
0360 auto subdetid = (id.subdetId());
0361 bool isPixel = subdetid == PixelSubdetector::PixelBarrel || subdetid == PixelSubdetector::PixelEndcap;
0362
0363 auto const& ltp = trajParams[h];
0364 auto ldir = ltp.momentum() / ltp.momentum().mag();
0365
0366 if (isPixel) {
0367
0368 const SiPixelRecHit* pixelRecHit = dynamic_cast<const SiPixelRecHit*>(recHit);
0369
0370 if (pixelRecHit != nullptr)
0371 processRec(*pixelRecHit, theClusterShape, ldir, *clusterShapeCache, hrpc);
0372 }
0373 }
0374 }
0375 }
0376
0377
0378 void PixelClusterShapeExtractor::analyze(edm::StreamID, const edm::Event& ev, const edm::EventSetup& es) const {
0379 if (hasSimHits) {
0380 LogTrace("MinBiasTracking") << " [ClusterShape] analyze simHits, recHits";
0381 analyzeSimHits(ev, es);
0382 }
0383
0384 if (hasRecTracks) {
0385 LogTrace("MinBiasTracking") << " [ClusterShape] analyze recHits on recTracks";
0386 analyzeRecTracks(ev, es);
0387 }
0388 }
0389
0390 DEFINE_FWK_MODULE(PixelClusterShapeExtractor);