PixelLowPtUtilities/bin/PixelClusterShapeExtractor.cc

0001 // VI January 2012: needs to be migrated to use cluster directly
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 // pixel
0043 #define exMax 10
0044 #define eyMax 15
0045
0046 namespace {
0047
0048   const std::set<unsigned> RelevantProcesses = {0};
0049   //const std::set<unsigned> RelevantProcesses = { 2, 7, 9, 11, 13, 15 };
0050
0051 }  // namespace
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   // Sim
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   // Rec
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   /// Tokens for ESconsumes
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;  // simulated pixel cluster
0112   std::vector<TH2F*> hrpc;  // reconstructed pixel cluster
0113 };
0114
0115 /*****************************************************************************/
0116 void PixelClusterShapeExtractor::init() {
0117   // Declare histograms
0118   char histName[256];
0119
0120   // pixel
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   // simulated
0171   for (auto h = hspc.begin(); h != hspc.end(); h++)
0172     (*h)->Write();
0173
0174   // reconstructed
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   // Outgoing?
0184   // very expensive....
0185   GlobalVector gvec = gdu.position() - GlobalPoint(0, 0, 0);
0186   LocalVector lvec = gdu.toLocal(gvec);
0187   LocalVector ldir = simHit.exitPoint() - simHit.entryPoint();
0188
0189   // cut on size as well (pixel is 285um thick...
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   // From a relevant process? primary or decay
0195   //bool isRelevant = (simHit.processType() == 2 ||
0196   //                   simHit.processType() == 4);
0197
0198   constexpr float ptCut2 = 0.2 * 0.2;  //  0.050*0.050;
0199   // Fast enough? pt > 50 MeV/c   FIXME (at least 200MeV....
0200   bool isFast = (simHit.momentumAtEntry().perp2() > ptCut2);
0201
0202   //std::cout << "isOutgoing = " << isOutgoing << ", isRelevant = " << simHit.processType() << ", isFast = " << isFast << std::endl;
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   //std::cout << "simHits.size() = " << simHits.size() << std::endl;
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     // Fill map
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     /* irrelevant
0294    auto const rh = *elem.second.rhit;
0295    auto const& topol = reinterpret_cast<const PixelGeomDetUnit*>(rh.detUnit())->specificTopology();
0296    auto const & cl = *rh.cluster();
0297    if (cl.minPixelCol()==0) continue;
0298    if (cl.maxPixelCol()+1==topol.ncolumns()) continue;
0299    if (cl.minPixelRow()==0) continue;
0300    if (cl.maxPixelRow()+1==topol.nrows()) continue;
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   // Get associator
0316   auto theHitAssociator = std::make_unique<TrackerHitAssociator>(ev, trackerHitAssociatorConfig_);
0317
0318   // Pixel hits
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   // Get tracks
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       // check that we are in the pixel
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         // Pixel
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);