Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-02-14 23:39:49 

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 "RecoPixelVertexing/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   size_t counter = 0, counter_2 = 0;
0296 
0297   for (SiPixelRecHitCollection::DataContainer::const_iterator recHit = recHits->begin(); recHit != recHits->end();
0298        recHit++)
0299     if (checkSimHits(*recHit, simHit, key)) {
0300       // Fill map
0301       if (simHitMap.count(key) == 0) {
0302         simHitMap[key] = &(*recHit);
0303       } else if (recHit->cluster()->size() > simHitMap[key]->cluster()->size())
0304         simHitMap[key] = &(*recHit);
0305       ++counter_2;
0306     }
0307 
0308   for (SiPixelRecHitCollection::DataContainer::const_iterator recHit = recHits->begin(); recHit != recHits->end();
0309        recHit++)
0310     if (checkSimHits(*recHit, simHit, key)) {
0311       // Check whether the present rechit is the largest
0312       if (&(*recHit) == simHitMap[key]) {
0313         processSim(*recHit, simHit, clusterShapeCache, hspc);
0314         ++counter;
0315       }
0316     }
0317   //    std::cout << "recHits->size() = " << recHits->size() << ", counter = " << counter
0318   //              << ", counter_2 = " << counter_2 << std::endl;
0319 }
0320 
0321 /*****************************************************************************/
0322 void ClusterShapeExtractor::processStripRecHits(const SiStripRecHit2DCollection::DataContainer* recHits) {
0323   map<pair<unsigned int, float>, const SiStripRecHit2D*> simHitMap;
0324 
0325   PSimHit simHit;
0326   pair<unsigned int, float> key;
0327 
0328   for (SiStripRecHit2DCollection::DataContainer::const_iterator recHit = recHits->begin(); recHit != recHits->end();
0329        recHit++)
0330     if (checkSimHits(*recHit, simHit, key)) {
0331       // Fill map
0332       if (simHitMap.count(key) == 0)
0333         simHitMap[key] = &(*recHit);
0334       else if (recHit->cluster()->amplitudes().size() > simHitMap[key]->cluster()->amplitudes().size())
0335         simHitMap[key] = &(*recHit);
0336     }
0337 
0338   for (SiStripRecHit2DCollection::DataContainer::const_iterator recHit = recHits->begin(); recHit != recHits->end();
0339        recHit++)
0340     if (checkSimHits(*recHit, simHit, key)) {
0341       // Check whether the present rechit is the largest
0342       if (&(*recHit) == simHitMap[key])
0343         processSim(*recHit, simHit, hssc);
0344     }
0345 }
0346 
0347 /*****************************************************************************/
0348 void ClusterShapeExtractor::processMatchedRecHits(const SiStripMatchedRecHit2DCollection::DataContainer* recHits) {
0349   map<pair<unsigned int, float>, const SiStripRecHit2D*> simHitMap;
0350   // VI very very quick fix
0351   std::vector<SiStripRecHit2D> cache;
0352   cache.reserve(2 * recHits->size());
0353   PSimHit simHit;
0354   pair<unsigned int, float> key;
0355 
0356   const SiStripRecHit2D* recHit;
0357 
0358   for (SiStripMatchedRecHit2DCollection::DataContainer::const_iterator matched = recHits->begin();
0359        matched != recHits->end();
0360        matched++) {
0361     cache.push_back(matched->monoHit());
0362     recHit = &cache.back();
0363     if (checkSimHits(*recHit, simHit, key)) {
0364       // Fill map
0365       if (simHitMap.count(key) == 0)
0366         simHitMap[key] = &(*recHit);
0367       else if (recHit->cluster()->amplitudes().size() > simHitMap[key]->cluster()->amplitudes().size())
0368         simHitMap[key] = &(*recHit);
0369     }
0370     cache.push_back(matched->stereoHit());
0371     recHit = &cache.back();
0372     if (checkSimHits(*recHit, simHit, key)) {
0373       // Fill map
0374       if (simHitMap.count(key) == 0)
0375         simHitMap[key] = &(*recHit);
0376       else if (recHit->cluster()->amplitudes().size() > simHitMap[key]->cluster()->amplitudes().size())
0377         simHitMap[key] = &(*recHit);
0378     }
0379   }
0380 
0381   for (SiStripMatchedRecHit2DCollection::DataContainer::const_iterator matched = recHits->begin();
0382        matched != recHits->end();
0383        matched++) {
0384     auto recHit = matched->monoHit();
0385     if (checkSimHits(recHit, simHit, key)) {
0386       // Check whether the present rechit is the largest
0387       if (recHit.omniCluster() == simHitMap[key]->omniCluster())
0388         processSim(recHit, simHit, hssc);
0389     }
0390 
0391     recHit = matched->stereoHit();
0392     if (checkSimHits(recHit, simHit, key)) {
0393       // Check whether the present rechit is the largest
0394       if (recHit.omniCluster() == simHitMap[key]->omniCluster())
0395         processSim(recHit, simHit, hssc);
0396     }
0397   }
0398 }
0399 
0400 /*****************************************************************************/
0401 void ClusterShapeExtractor::analyzeSimHits(const edm::Event& ev, const edm::EventSetup& es) {
0402   // Get associator
0403   theHitAssociator.reset(new TrackerHitAssociator(ev, trackerHitAssociatorConfig_));
0404 
0405   // Pixel hits
0406   {
0407     edm::Handle<SiPixelRecHitCollection> coll;
0408     //ev.getByLabel("siPixelRecHits", coll);
0409     ev.getByToken(pixelRecHits_token, coll);
0410 
0411     edm::Handle<SiPixelClusterShapeCache> clusterShapeCache;
0412     ev.getByToken(clusterShapeCache_token, clusterShapeCache);
0413 
0414     const SiPixelRecHitCollection::DataContainer* recHits = &coll.product()->data();
0415     processPixelRecHits(recHits, *clusterShapeCache);
0416   }
0417 
0418   /*
0419   // Strip hits
0420   { // rphi and stereo
0421     vector<edm::Handle<SiStripRecHit2DCollection> > colls;
0422     ev.getManyByType(colls);
0423   
0424     for(vector<edm::Handle<SiStripRecHit2DCollection> >::const_iterator
0425           coll = colls.begin(); coll!= colls.end(); coll++)
0426     {
0427       const SiStripRecHit2DCollection::DataContainer * recHits =
0428             & (*coll).product()->data();
0429       processStripRecHits(recHits);
0430     }
0431   }
0432 
0433   // Matched strip hits
0434   { // matched
0435   vector<edm::Handle<SiStripMatchedRecHit2DCollection> > colls;
0436   ev.getManyByType(colls);
0437 
0438   for(vector<edm::Handle<SiStripMatchedRecHit2DCollection> >::const_iterator
0439         coll = colls.begin(); coll!= colls.end(); coll++)
0440   {
0441     const SiStripMatchedRecHit2DCollection::DataContainer * recHits =
0442           & (*coll).product()->data();
0443     processMatchedRecHits(recHits);
0444   }
0445   }
0446   */
0447 }
0448 
0449 /*****************************************************************************/
0450 void ClusterShapeExtractor::analyzeRecTracks(const edm::Event& ev, const edm::EventSetup& es) {
0451   // Get trajectory collection
0452   edm::Handle<vector<Trajectory>> trajeHandle;
0453   ev.getByLabel(trackProducer, trajeHandle);
0454   const vector<Trajectory>& trajeCollection = *(trajeHandle.product());
0455 
0456   edm::Handle<SiPixelClusterShapeCache> clusterShapeCache;
0457   ev.getByToken(clusterShapeCache_token, clusterShapeCache);
0458 
0459   // Take all trajectories
0460   for (vector<Trajectory>::const_iterator trajectory = trajeCollection.begin(); trajectory != trajeCollection.end();
0461        trajectory++)
0462     for (vector<TrajectoryMeasurement>::const_iterator meas = trajectory->measurements().begin();
0463          meas != trajectory->measurements().end();
0464          meas++) {
0465       const TrackingRecHit* recHit = meas->recHit()->hit();
0466       DetId id = recHit->geographicalId();
0467 
0468       if (recHit->isValid()) {
0469         LocalVector ldir = meas->updatedState().localDirection();
0470 
0471         if (GeomDetEnumerators::isTrackerPixel(theTracker->geomDetSubDetector(id.subdetId()))) {
0472           // Pixel
0473           const SiPixelRecHit* pixelRecHit = dynamic_cast<const SiPixelRecHit*>(recHit);
0474 
0475           if (pixelRecHit != 0)
0476             processRec(*pixelRecHit, ldir, *clusterShapeCache, hrpc);
0477         } else if (GeomDetEnumerators::isTrackerStrip(theTracker->geomDetSubDetector(id.subdetId()))) {
0478           // Strip
0479           const SiStripMatchedRecHit2D* stripMatchedRecHit = dynamic_cast<const SiStripMatchedRecHit2D*>(recHit);
0480           const ProjectedSiStripRecHit2D* stripProjectedRecHit = dynamic_cast<const ProjectedSiStripRecHit2D*>(recHit);
0481           const SiStripRecHit2D* stripRecHit = dynamic_cast<const SiStripRecHit2D*>(recHit);
0482 
0483           if (stripMatchedRecHit != 0) {
0484             processRec(stripMatchedRecHit->monoHit(), ldir, hrsc);
0485             processRec(stripMatchedRecHit->stereoHit(), ldir, hrsc);
0486           }
0487 
0488           if (stripProjectedRecHit != 0)
0489             processRec(stripProjectedRecHit->originalHit(), ldir, hrsc);
0490 
0491           if (stripRecHit != 0)
0492             processRec(*stripRecHit, ldir, hrsc);
0493         }
0494       }
0495     }
0496 }
0497 
0498 /*****************************************************************************/
0499 void ClusterShapeExtractor::analyze(const edm::Event& ev, const edm::EventSetup& es) {
0500   if (hasSimHits) {
0501     LogTrace("MinBiasTracking") << " [ClusterShape] analyze simHits, recHits";
0502     analyzeSimHits(ev, es);
0503   }
0504 
0505   if (hasRecTracks) {
0506     LogTrace("MinBiasTracking") << " [ClusterShape] analyze recHits on recTracks";
0507     analyzeRecTracks(ev, es);
0508   }
0509 }
0510 
0511 DEFINE_FWK_MODULE(ClusterShapeExtractor);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!