PixelLowPtUtilities/test/ClusterShapeExtractor.cc

0001 // VI January 2012: needs to be migrated to use cluster directly
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 // pixel
0042 #define exMax 10
0043 #define eyMax 15
0044
0045 // strip
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   // Sim
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   // Rec
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   // member data
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;  // simulated pixel cluster
0104   vector<TH1F*> hssc;  // simulated strip cluster
0105
0106   vector<TH2F*> hrpc;  // reconstructed pixel cluster
0107   vector<TH1F*> hrsc;  // reconstructed strip cluster
0108 };
0109
0110 /*****************************************************************************/
0111 void ClusterShapeExtractor::beginRun(const edm::Run& run, const edm::EventSetup& es) {
0112   // Get tracker geometry
0113   theTracker = &es.getData(geomToken);
0114
0115   // Get the cluster shape
0116   theClusterShape = &es.getData(clusterShapeToken);
0117
0118   // Declare histograms
0119   char histName[256];
0120
0121   // pixel
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   // strip
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   // simulated
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   // reconstructed
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   //  delete theClusterShape;
0198 }
0199
0200 /*****************************************************************************/
0201 bool ClusterShapeExtractor::isSuitable(const PSimHit& simHit) {
0202   // Outgoing?
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   //static const std::set<unsigned> RelevantProcesses = { 2, 7, 9, 11, 13, 15 };
0215   const bool isRelevant = RelevantProcesses.count(simHit.processType());
0216   // From a relevant process? primary or decay
0217   //bool isRelevant = (simHit.processType() == 2 ||
0218   //                   simHit.processType() == 4);
0219
0220   // Fast enough? pt > 50 MeV/c
0221   bool isFast = (simHit.momentumAtEntry().perp() > 0.050);
0222
0223   //std::cout << "isOutgoing = " << isOutgoing << ", isRelevant = " << simHit.processType() << ", isFast = " << isFast << std::endl;
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   //std::cout << "simHits.size() = " << simHits.size() << std::endl;
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       // Fill map
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       // Check whether the present rechit is the largest
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       // Fill map
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       // Check whether the present rechit is the largest
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   // VI very very quick fix
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       // Fill map
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       // Fill map
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       // Check whether the present rechit is the largest
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       // Check whether the present rechit is the largest
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   // Get associator
0398   theHitAssociator.reset(new TrackerHitAssociator(ev, trackerHitAssociatorConfig_));
0399
0400   // Pixel hits
0401   {
0402     edm::Handle<SiPixelRecHitCollection> coll;
0403     //ev.getByLabel("siPixelRecHits", coll);
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   // Strip hits
0415   { // rphi and stereo
0416     vector<edm::Handle<SiStripRecHit2DCollection> > colls;
0417     ev.getManyByType(colls);
0418
0419     for(vector<edm::Handle<SiStripRecHit2DCollection> >::const_iterator
0420           coll = colls.begin(); coll!= colls.end(); coll++)
0421     {
0422       const SiStripRecHit2DCollection::DataContainer * recHits =
0423             & (*coll).product()->data();
0424       processStripRecHits(recHits);
0425     }
0426   }
0427
0428   // Matched strip hits
0429   { // matched
0430   vector<edm::Handle<SiStripMatchedRecHit2DCollection> > colls;
0431   ev.getManyByType(colls);
0432
0433   for(vector<edm::Handle<SiStripMatchedRecHit2DCollection> >::const_iterator
0434         coll = colls.begin(); coll!= colls.end(); coll++)
0435   {
0436     const SiStripMatchedRecHit2DCollection::DataContainer * recHits =
0437           & (*coll).product()->data();
0438     processMatchedRecHits(recHits);
0439   }
0440   }
0441   */
0442 }
0443
0444 /*****************************************************************************/
0445 void ClusterShapeExtractor::analyzeRecTracks(const edm::Event& ev, const edm::EventSetup& es) {
0446   // Get trajectory collection
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   // Take all trajectories
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           // Pixel
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           // Strip
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);