File indexing completed on 2024-04-06 12:26:16
0001
0002 #include <map>
0003 #include <vector>
0004 #include <algorithm>
0005
0006
0007 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0008 #include "CommonTools/Utils/interface/TFileDirectory.h"
0009 #include "DataFormats/Common/interface/DetSetVector.h"
0010 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "DataFormats/DetId/interface/DetId.h"
0013 #include "DataFormats/Phase2TrackerCluster/interface/Phase2TrackerCluster1D.h"
0014 #include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h"
0015 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0016 #include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
0017 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0018 #include "FWCore/Framework/interface/ConsumesCollector.h"
0019 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0020 #include "FWCore/Framework/interface/ESHandle.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/Frameworkfwd.h"
0023 #include "FWCore/Framework/interface/MakerMacros.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/ServiceRegistry/interface/Service.h"
0026 #include "FWCore/ServiceRegistry/interface/Service.h"
0027 #include "FWCore/Utilities/interface/InputTag.h"
0028 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0029 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0030 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0031 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0032 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0033 #include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h"
0034 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0035 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0036
0037
0038 #include <TH2F.h>
0039 #include <TH1F.h>
0040 #include <THStack.h>
0041
0042 struct ClusterHistos {
0043
0044
0045
0046 TH1D* numberClusters[3];
0047 TH1D* clusterSize[3];
0048
0049 TH2F* globalPosXY[3];
0050 TH2F* localPosXY[3];
0051
0052 TH1F* deltaX[3];
0053 TH1F* deltaY[3];
0054 TH1F* deltaX_P[3];
0055 TH1F* deltaY_P[3];
0056
0057 TH1D* primarySimHits[3];
0058 TH1D* otherSimHits[3];
0059 };
0060
0061 class Phase2TrackerClusterizerValidation : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0062 public:
0063 typedef std::map<unsigned int, std::vector<PSimHit> > SimHitsMap;
0064 typedef std::map<unsigned int, SimTrack> SimTracksMap;
0065
0066 explicit Phase2TrackerClusterizerValidation(const edm::ParameterSet&);
0067 ~Phase2TrackerClusterizerValidation();
0068 void beginJob() override;
0069 void analyze(const edm::Event&, const edm::EventSetup&) override;
0070
0071 private:
0072 std::map<unsigned int, ClusterHistos>::iterator createLayerHistograms(unsigned int);
0073 std::vector<unsigned int> getSimTrackId(const edm::Handle<edm::DetSetVector<PixelDigiSimLink> >&,
0074 const DetId&,
0075 unsigned int);
0076
0077 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> esTokenGeom_;
0078 const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> esTokenTopo_;
0079 const edm::EDGetTokenT<Phase2TrackerCluster1DCollectionNew> tokenClusters_;
0080 const edm::EDGetTokenT<edm::DetSetVector<PixelDigiSimLink> > tokenLinks_;
0081 const edm::EDGetTokenT<edm::PSimHitContainer> tokenSimHitsB_;
0082 const edm::EDGetTokenT<edm::PSimHitContainer> tokenSimHitsE_;
0083 const edm::EDGetTokenT<edm::SimTrackContainer> tokenSimTracks_;
0084
0085 const bool catECasRings_;
0086 const double simtrackminpt_;
0087
0088 TH2F* trackerLayout_;
0089 TH2F* trackerLayoutXY_;
0090 TH2F* trackerLayoutXYBar_;
0091 TH2F* trackerLayoutXYEC_;
0092
0093 std::map<unsigned int, ClusterHistos> histograms_;
0094 };
0095
0096 Phase2TrackerClusterizerValidation::Phase2TrackerClusterizerValidation(const edm::ParameterSet& conf)
0097 : esTokenGeom_(esConsumes()),
0098 esTokenTopo_(esConsumes()),
0099 tokenClusters_(consumes<Phase2TrackerCluster1DCollectionNew>(conf.getParameter<edm::InputTag>("src"))),
0100 tokenLinks_(consumes<edm::DetSetVector<PixelDigiSimLink> >(conf.getParameter<edm::InputTag>("links"))),
0101 tokenSimHitsB_(consumes<edm::PSimHitContainer>(conf.getParameter<edm::InputTag>("simhitsbarrel"))),
0102 tokenSimHitsE_(consumes<edm::PSimHitContainer>(conf.getParameter<edm::InputTag>("simhitsendcap"))),
0103 tokenSimTracks_(consumes<edm::SimTrackContainer>(conf.getParameter<edm::InputTag>("simtracks"))),
0104 catECasRings_(conf.getParameter<bool>("ECasRings")),
0105 simtrackminpt_(conf.getParameter<double>("SimTrackMinPt")) {
0106 usesResource(TFileService::kSharedResource);
0107 }
0108
0109 Phase2TrackerClusterizerValidation::~Phase2TrackerClusterizerValidation() = default;
0110
0111 void Phase2TrackerClusterizerValidation::beginJob() {
0112 edm::Service<TFileService> fs;
0113 fs->file().cd("/");
0114 TFileDirectory td = fs->mkdir("Common");
0115
0116 trackerLayout_ = td.make<TH2F>("RVsZ", "R vs. z position", 6000, -300.0, 300.0, 1200, 0.0, 120.0);
0117 trackerLayoutXY_ = td.make<TH2F>("XVsY", "x vs. y position", 2400, -120.0, 120.0, 2400, -120.0, 120.0);
0118 trackerLayoutXYBar_ = td.make<TH2F>("XVsYBar", "x vs. y position", 2400, -120.0, 120.0, 2400, -120.0, 120.0);
0119 trackerLayoutXYEC_ = td.make<TH2F>("XVsYEC", "x vs. y position", 2400, -120.0, 120.0, 2400, -120.0, 120.0);
0120 }
0121
0122 void Phase2TrackerClusterizerValidation::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0123
0124
0125
0126
0127
0128 edm::Handle<Phase2TrackerCluster1DCollectionNew> clusters;
0129 event.getByToken(tokenClusters_, clusters);
0130
0131
0132 edm::Handle<edm::DetSetVector<PixelDigiSimLink> > pixelSimLinks;
0133 event.getByToken(tokenLinks_, pixelSimLinks);
0134
0135
0136 edm::Handle<edm::PSimHitContainer> simHitsRaw[2];
0137 event.getByToken(tokenSimHitsB_, simHitsRaw[0]);
0138 event.getByToken(tokenSimHitsE_, simHitsRaw[1]);
0139
0140
0141 edm::Handle<edm::SimTrackContainer> simTracksRaw;
0142 event.getByToken(tokenSimTracks_, simTracksRaw);
0143
0144
0145 const TrackerGeometry* tkGeom = &eventSetup.getData(esTokenGeom_);
0146 const TrackerTopology* tTopo = &eventSetup.getData(esTokenTopo_);
0147
0148
0149
0150
0151
0152
0153 SimTracksMap simTracks;
0154 for (edm::SimTrackContainer::const_iterator simTrackIt(simTracksRaw->begin()); simTrackIt != simTracksRaw->end();
0155 ++simTrackIt) {
0156 if (simTrackIt->momentum().pt() > simtrackminpt_) {
0157 simTracks.emplace(simTrackIt->trackId(), *simTrackIt);
0158 }
0159 }
0160
0161
0162
0163
0164
0165
0166 std::map<unsigned int, unsigned int> nClusters[3];
0167 std::map<unsigned int, unsigned int> nPrimarySimHits[3];
0168 std::map<unsigned int, unsigned int> nOtherSimHits[3];
0169
0170
0171 for (Phase2TrackerCluster1DCollectionNew::const_iterator DSViter = clusters->begin(); DSViter != clusters->end();
0172 ++DSViter) {
0173
0174 unsigned int rawid(DSViter->detId());
0175 DetId detId(rawid);
0176 unsigned int layer = (tTopo->side(detId) != 0) * 1000;
0177 if (!layer) {
0178 layer += tTopo->layer(detId);
0179 } else {
0180 layer += (catECasRings_ ? tTopo->tidRing(detId) * 10 : tTopo->layer(detId));
0181 }
0182 TrackerGeometry::ModuleType mType = tkGeom->getDetectorType(detId);
0183 unsigned int det = 0;
0184 if (mType == TrackerGeometry::ModuleType::Ph2PSP) {
0185 det = 1;
0186 } else if (mType == TrackerGeometry::ModuleType::Ph2PSS || mType == TrackerGeometry::ModuleType::Ph2SS) {
0187 det = 2;
0188 } else {
0189 std::cout << "UNKNOWN DETECTOR TYPE!" << std::endl;
0190 }
0191
0192
0193 const GeomDetUnit* geomDetUnit(tkGeom->idToDetUnit(detId));
0194 if (!geomDetUnit)
0195 continue;
0196
0197
0198 auto nhitit(nClusters[det].find(layer));
0199 if (nhitit == nClusters[det].end()) {
0200 nClusters[det].emplace(layer, 0);
0201 nPrimarySimHits[det].emplace(layer, 0);
0202 nOtherSimHits[det].emplace(layer, 0);
0203 }
0204
0205
0206 std::map<unsigned int, ClusterHistos>::iterator histogramLayer(histograms_.find(layer));
0207 if (histogramLayer == histograms_.end())
0208 histogramLayer = createLayerHistograms(layer);
0209
0210
0211 for (edmNew::DetSet<Phase2TrackerCluster1D>::const_iterator clustIt = DSViter->begin(); clustIt != DSViter->end();
0212 ++clustIt) {
0213
0214 MeasurementPoint mpClu(clustIt->center(), clustIt->column() + 0.5);
0215 Local3DPoint localPosClu = geomDetUnit->topology().localPosition(mpClu);
0216 Global3DPoint globalPosClu = geomDetUnit->surface().toGlobal(localPosClu);
0217
0218
0219 std::vector<unsigned int> clusterSimTrackIds;
0220 for (unsigned int i(0); i < clustIt->size(); ++i) {
0221 unsigned int channel(Phase2TrackerDigi::pixelToChannel(clustIt->firstRow() + i, clustIt->column()));
0222 std::vector<unsigned int> simTrackIds(getSimTrackId(pixelSimLinks, detId, channel));
0223 for (auto it : simTrackIds) {
0224 bool add = true;
0225 for (unsigned int j = 0; j < clusterSimTrackIds.size(); ++j) {
0226
0227 if (it == clusterSimTrackIds.at(j))
0228 add = false;
0229 }
0230 if (add)
0231 clusterSimTrackIds.push_back(it);
0232 }
0233 }
0234 std::sort(clusterSimTrackIds.begin(), clusterSimTrackIds.end());
0235
0236
0237
0238
0239 const PSimHit* simhit = 0;
0240 float minx = 10000;
0241 for (unsigned int simhitidx = 0; simhitidx < 2; ++simhitidx) {
0242 for (auto simhitIt : *simHitsRaw[simhitidx]) {
0243 if (rawid == simhitIt.detUnitId()) {
0244
0245 auto it = std::lower_bound(clusterSimTrackIds.begin(), clusterSimTrackIds.end(), simhitIt.trackId());
0246 if (it != clusterSimTrackIds.end() && *it == simhitIt.trackId()) {
0247 if (!simhit || fabs(simhitIt.localPosition().x() - localPosClu.x()) < minx) {
0248 minx = fabs(simhitIt.localPosition().x() - localPosClu.x());
0249 simhit = &simhitIt;
0250 }
0251 }
0252 }
0253 }
0254 }
0255 if (!simhit)
0256 continue;
0257
0258
0259 auto simTrackIt(simTracks.find(simhit->trackId()));
0260 if (simTrackIt == simTracks.end())
0261 continue;
0262
0263
0264
0265
0266
0267
0268 ++(nClusters[det].at(layer));
0269 ++(nOtherSimHits[det].at(layer));
0270
0271
0272 histogramLayer->second.clusterSize[det]->Fill(clustIt->size());
0273
0274
0275 trackerLayout_->Fill(globalPosClu.z(), globalPosClu.perp());
0276 trackerLayoutXY_->Fill(globalPosClu.x(), globalPosClu.y());
0277 if (layer < 1000) {
0278 trackerLayoutXYBar_->Fill(globalPosClu.x(), globalPosClu.y());
0279 } else {
0280 trackerLayoutXYEC_->Fill(globalPosClu.x(), globalPosClu.y());
0281 }
0282
0283 histogramLayer->second.localPosXY[det]->Fill(localPosClu.x(), localPosClu.y());
0284 if (layer < 1000) {
0285 histogramLayer->second.globalPosXY[det]->Fill(globalPosClu.z(), globalPosClu.perp());
0286 } else {
0287 histogramLayer->second.globalPosXY[det]->Fill(globalPosClu.x(), globalPosClu.y());
0288 }
0289
0290
0291 Local3DPoint localPosHit(simhit->localPosition());
0292
0293 histogramLayer->second.deltaX[det]->Fill(localPosClu.x() - localPosHit.x());
0294 histogramLayer->second.deltaY[det]->Fill(localPosClu.y() - localPosHit.y());
0295
0296
0297 unsigned int procT(simhit->processType());
0298 if (simTrackIt->second.vertIndex() == 0 and
0299 (procT == 2 || procT == 7 || procT == 9 || procT == 11 || procT == 13 || procT == 15)) {
0300 ++(nPrimarySimHits[det].at(layer));
0301 --(nOtherSimHits[det].at(layer));
0302 histogramLayer->second.deltaX_P[det]->Fill(localPosClu.x() - localPosHit.x());
0303 histogramLayer->second.deltaY_P[det]->Fill(localPosClu.y() - localPosHit.y());
0304 }
0305 }
0306 }
0307
0308
0309 for (unsigned int det = 1; det < 3; ++det) {
0310 for (auto it : nClusters[det]) {
0311 auto histogramLayer(histograms_.find(it.first));
0312 if (histogramLayer == histograms_.end())
0313 std::cout << "*** SL *** No histogram for an existing counter! This should not happen!" << std::endl;
0314 histogramLayer->second.numberClusters[det]->Fill(it.second);
0315 }
0316 for (auto it : nPrimarySimHits[det]) {
0317 auto histogramLayer(histograms_.find(it.first));
0318 if (histogramLayer == histograms_.end())
0319 std::cout << "*** SL *** No histogram for an existing counter! This should not happen!" << std::endl;
0320 histogramLayer->second.primarySimHits[det]->Fill(it.second);
0321 }
0322 for (auto it : nOtherSimHits[det]) {
0323 auto histogramLayer(histograms_.find(it.first));
0324 if (histogramLayer == histograms_.end())
0325 std::cout << "*** SL *** No histogram for an existing counter! This should not happen!" << std::endl;
0326 histogramLayer->second.otherSimHits[det]->Fill(it.second);
0327 }
0328 }
0329 }
0330
0331
0332 std::map<unsigned int, ClusterHistos>::iterator Phase2TrackerClusterizerValidation::createLayerHistograms(
0333 unsigned int ival) {
0334 std::ostringstream fname1, fname2;
0335
0336 edm::Service<TFileService> fs;
0337 fs->file().cd("/");
0338
0339 std::string tag;
0340 unsigned int id;
0341 if (ival < 1000) {
0342 id = ival;
0343 fname1 << "Barrel";
0344 fname2 << "Layer_" << id;
0345 tag = "_layer_";
0346 } else {
0347 int side = ival / 1000;
0348 id = ival - side * 1000;
0349 if (ival > 10) {
0350 id /= 10;
0351
0352 fname1 << "EndCap";
0353 fname2 << "Ring_" << id;
0354 tag = "_ring_";
0355 } else {
0356 id = ival;
0357
0358 fname1 << "EndCap";
0359 fname2 << "Disk_" << id;
0360 tag = "_disk_";
0361 }
0362 }
0363
0364 TFileDirectory td1 = fs->mkdir(fname1.str().c_str());
0365 TFileDirectory td = td1.mkdir(fname2.str().c_str());
0366
0367 ClusterHistos local_histos;
0368
0369 std::ostringstream histoName;
0370
0371
0372
0373
0374
0375 histoName.str("");
0376 histoName << "Number_Clusters_Pixel" << tag.c_str() << id;
0377 local_histos.numberClusters[1] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 50, 0., 0.);
0378 histoName.str("");
0379 histoName << "Number_Clusters_Strip" << tag.c_str() << id;
0380 local_histos.numberClusters[2] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 50, 0., 0.);
0381
0382
0383
0384
0385
0386 histoName.str("");
0387 histoName << "Cluster_Size_Pixel" << tag.c_str() << id;
0388 local_histos.clusterSize[1] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 21, -0.5, 20.5);
0389 histoName.str("");
0390 histoName << "Cluster_Size_Strip" << tag.c_str() << id;
0391 local_histos.clusterSize[2] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 21, -0.5, 20.5);
0392
0393
0394
0395
0396
0397 histoName.str("");
0398 histoName << "Local_Position_XY_Pixel" << tag.c_str() << id;
0399 local_histos.localPosXY[1] =
0400 td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 500, 0., 0., 500, 0., 0.);
0401
0402 histoName.str("");
0403 histoName << "Local_Position_XY_Strip" << tag.c_str() << id;
0404 local_histos.localPosXY[2] =
0405 td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 500, 0., 0., 500, 0., 0.);
0406
0407 histoName.str("");
0408 histoName << "Global_Position_XY_Pixel" << tag.c_str() << id;
0409 local_histos.globalPosXY[1] =
0410 td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 500, -120.0, 120.0, 2400, -120.0, 120.0);
0411
0412 histoName.str("");
0413 histoName << "Global_Position_XY_Strip" << tag.c_str() << id;
0414 local_histos.globalPosXY[2] =
0415 td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 500, -120.0, 120.0, 2400, -120.0, 120.0);
0416
0417
0418
0419
0420
0421 histoName.str("");
0422 histoName << "Delta_X_Pixel" << tag.c_str() << id;
0423 local_histos.deltaX[1] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -0.02, 0.02);
0424
0425 histoName.str("");
0426 histoName << "Delta_X_Strip" << tag.c_str() << id;
0427 local_histos.deltaX[2] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -0.02, 0.02);
0428
0429 histoName.str("");
0430 histoName << "Delta_Y_Pixel" << tag.c_str() << id;
0431 local_histos.deltaY[1] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -0.2, 0.2);
0432
0433 histoName.str("");
0434 histoName << "Delta_Y_Strip" << tag.c_str() << id;
0435 local_histos.deltaY[2] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -3., 3.);
0436
0437
0438
0439
0440
0441 histoName.str("");
0442 histoName << "Delta_X_P" << tag.c_str() << id;
0443 local_histos.deltaX_P[1] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -0.02, 0.02);
0444
0445 histoName.str("");
0446 histoName << "Delta_X_P" << tag.c_str() << id;
0447 local_histos.deltaX_P[2] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -0.02, 0.02);
0448
0449 histoName.str("");
0450 histoName << "Delta_Y_P" << tag.c_str() << id;
0451 local_histos.deltaY_P[1] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -0.2, 0.2);
0452
0453 histoName.str("");
0454 histoName << "Delta_Y_P" << tag.c_str() << id;
0455 local_histos.deltaY_P[2] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -3., 3.);
0456
0457
0458
0459
0460
0461 histoName.str("");
0462 histoName << "Primary_Digis_Pixel" << tag.c_str() << id;
0463 local_histos.primarySimHits[1] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
0464 histoName.str("");
0465 histoName << "Primary_Digis_Strip" << tag.c_str() << id;
0466 local_histos.primarySimHits[2] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
0467
0468 histoName.str("");
0469 histoName << "Other_Digis_Pixel" << tag.c_str() << id;
0470 local_histos.otherSimHits[1] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
0471 histoName.str("");
0472 histoName << "Other_Digis_Strip" << tag.c_str() << id;
0473 local_histos.otherSimHits[2] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
0474
0475
0476
0477
0478
0479 std::pair<std::map<unsigned int, ClusterHistos>::iterator, bool> insertedIt(
0480 histograms_.insert(std::make_pair(ival, local_histos)));
0481 fs->file().cd("/");
0482
0483 return insertedIt.first;
0484 }
0485
0486 std::vector<unsigned int> Phase2TrackerClusterizerValidation::getSimTrackId(
0487 const edm::Handle<edm::DetSetVector<PixelDigiSimLink> >& pixelSimLinks, const DetId& detId, unsigned int channel) {
0488 std::vector<unsigned int> retvec;
0489 edm::DetSetVector<PixelDigiSimLink>::const_iterator DSViter(pixelSimLinks->find(detId));
0490 if (DSViter == pixelSimLinks->end())
0491 return retvec;
0492 for (edm::DetSet<PixelDigiSimLink>::const_iterator it = DSViter->data.begin(); it != DSViter->data.end(); ++it) {
0493 if (channel == it->channel()) {
0494 retvec.push_back(it->SimTrackId());
0495 }
0496 }
0497 return retvec;
0498 }
0499
0500 DEFINE_FWK_MODULE(Phase2TrackerClusterizerValidation);