Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:38:02 

0001 // -*- C++ -*-
0002 //
0003 // Package:    trackJet/DeepCoreSeedGenerator
0004 // Class:      DeepCoreSeedGenerator
0005 //
0006 /**\class DeepCoreSeedGenerator DeepCoreSeedGenerator.cc trackJet/DeepCoreSeedGenerator/plugins/DeepCoreSeedGenerator.cc
0007 
0008  Description: [one line class summary]
0009 
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Valerio Bertacchi
0015 //         Created:  Mon, 18 Dec 2017 16:35:04 GMT
0016 //
0017 //
0018 
0019 // system include files
0020 
0021 #include <memory>
0022 
0023 // user include files
0024 
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/stream/EDProducer.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/EventSetup.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030 #include "FWCore/Framework/interface/ESHandle.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 
0033 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0034 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0035 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0036 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0037 
0038 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0039 #include "DataFormats/Common/interface/Handle.h"
0040 #include "DataFormats/Common/interface/DetSetVector.h"
0041 #include "DataFormats/Common/interface/DetSet.h"
0042 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0043 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0044 #include "DataFormats/VertexReco/interface/Vertex.h"
0045 #include "DataFormats/TrackReco/interface/Track.h"
0046 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0047 #include "DataFormats/JetReco/interface/Jet.h"
0048 #include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
0049 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0050 #include "DataFormats/Math/interface/Point3D.h"
0051 #include "DataFormats/Math/interface/Vector3D.h"
0052 #include "DataFormats/Candidate/interface/Candidate.h"
0053 
0054 #include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h"
0055 #include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h"
0056 
0057 #include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h"
0058 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0059 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0060 
0061 #include "MagneticField/Engine/interface/MagneticField.h"
0062 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0063 
0064 #include "FWCore/ServiceRegistry/interface/Service.h"
0065 
0066 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0067 
0068 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0069 
0070 #include "PhysicsTools/TensorFlow/interface/TensorFlow.h"
0071 
0072 //
0073 // class declaration
0074 //
0075 
0076 class DeepCoreSeedGenerator : public edm::stream::EDProducer<edm::GlobalCache<tensorflow::SessionCache>> {
0077 public:
0078   explicit DeepCoreSeedGenerator(const edm::ParameterSet&, const tensorflow::SessionCache*);
0079   ~DeepCoreSeedGenerator() override;
0080 
0081   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0082   // A pointer to a cluster and a list of tracks on it
0083 
0084   // static methods for handling the global cache
0085   static std::unique_ptr<tensorflow::SessionCache> initializeGlobalCache(const edm::ParameterSet&);
0086   static void globalEndJob(tensorflow::SessionCache*);
0087 
0088   double jetPt_;
0089   double jetEta_;
0090   double pitchX_ = 0.01;              //100 um (pixel pitch in X)
0091   double pitchY_ = 0.015;             //150 um (pixel pitch in Y)
0092   static constexpr int jetDimX = 30;  //pixel dimension of NN window on layer2
0093   static constexpr int jetDimY = 30;  //pixel dimension of NN window on layer2
0094   static constexpr int Nlayer = 4;    //Number of layer used in DeepCore
0095   static constexpr int Nover = 3;     //Max number of tracks recorded per pixel
0096   static constexpr int Npar = 5;      //Number of track parameter
0097 
0098   // DeepCore 2.2.1 thresholds delta
0099   double dth[3] = {0.0, 0.15, 0.3};
0100   double dthl[3] = {0.1, 0.05, 0};
0101 
0102 private:
0103   void produce(edm::Event&, const edm::EventSetup&) override;
0104 
0105   // ----------member data ---------------------------
0106   std::string propagatorName_;
0107   const GlobalTrackingGeometry* geometry_ = nullptr;
0108 
0109   edm::EDGetTokenT<std::vector<reco::Vertex>> vertices_;
0110   edm::EDGetTokenT<edmNew::DetSetVector<SiPixelCluster>> pixelClusters_;
0111   edm::Handle<edmNew::DetSetVector<SiPixelCluster>> inputPixelClusters_;
0112   edm::EDGetTokenT<edm::View<reco::Candidate>> cores_;
0113   const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> geometryToken_;
0114   const edm::ESGetToken<PixelClusterParameterEstimator, TkPixelCPERecord> pixelCPEToken_;
0115   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0116 
0117   double ptMin_;
0118   double deltaR_;
0119   double chargeFracMin_;
0120   double centralMIPCharge_;
0121   std::string weightfilename_;
0122   std::vector<std::string> inputTensorName_;
0123   std::vector<std::string> outputTensorName_;
0124   double probThr_;
0125   const tensorflow::Session* session_;
0126 
0127   std::pair<bool, Basic3DVector<float>> findIntersection(const GlobalVector&,
0128                                                          const reco::Candidate::Point&,
0129                                                          const GeomDet*);
0130 
0131   void fillPixelMatrix(const SiPixelCluster&,
0132                        int,
0133                        Point3DBase<float, LocalTag>,
0134                        const GeomDet*,
0135                        tensorflow::NamedTensorList);  //if not working,: args=2 auto
0136 
0137   std::pair<int, int> local2Pixel(double, double, const GeomDet*);
0138 
0139   LocalPoint pixel2Local(int, int, const GeomDet*);
0140 
0141   int pixelFlipper(const GeomDet*);
0142 
0143   const GeomDet* DetectorSelector(int,
0144                                   const reco::Candidate&,
0145                                   GlobalVector,
0146                                   const reco::Vertex&,
0147                                   const TrackerTopology* const,
0148                                   const edmNew::DetSetVector<SiPixelCluster>&);
0149 
0150   std::vector<GlobalVector> splittedClusterDirections(
0151       const reco::Candidate&,
0152       const TrackerTopology* const,
0153       const PixelClusterParameterEstimator*,
0154       const reco::Vertex&,
0155       int,
0156       const edmNew::DetSetVector<SiPixelCluster>&);  //if not working,: args=2 auto
0157 
0158   std::pair<double[jetDimX][jetDimY][Nover][Npar], double[jetDimX][jetDimY][Nover]> SeedEvaluation(
0159       tensorflow::NamedTensorList, std::vector<std::string>);
0160 };
0161 
0162 DeepCoreSeedGenerator::DeepCoreSeedGenerator(const edm::ParameterSet& iConfig, const tensorflow::SessionCache* cache)
0163     : vertices_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0164       pixelClusters_(
0165           consumes<edmNew::DetSetVector<SiPixelCluster>>(iConfig.getParameter<edm::InputTag>("pixelClusters"))),
0166       cores_(consumes<edm::View<reco::Candidate>>(iConfig.getParameter<edm::InputTag>("cores"))),
0167       geometryToken_(esConsumes()),
0168       pixelCPEToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("pixelCPE")))),
0169       topoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>()),
0170       ptMin_(iConfig.getParameter<double>("ptMin")),
0171       deltaR_(iConfig.getParameter<double>("deltaR")),
0172       chargeFracMin_(iConfig.getParameter<double>("chargeFractionMin")),
0173       centralMIPCharge_(iConfig.getParameter<double>("centralMIPCharge")),
0174       weightfilename_(iConfig.getParameter<edm::FileInPath>("weightFile").fullPath()),
0175       inputTensorName_(iConfig.getParameter<std::vector<std::string>>("inputTensorName")),
0176       outputTensorName_(iConfig.getParameter<std::vector<std::string>>("outputTensorName")),
0177       probThr_(iConfig.getParameter<double>("probThr")),
0178       session_(cache->getSession())
0179 
0180 {
0181   produces<TrajectorySeedCollection>();
0182   produces<reco::TrackCollection>();
0183 }
0184 
0185 DeepCoreSeedGenerator::~DeepCoreSeedGenerator() {}
0186 
0187 void DeepCoreSeedGenerator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0188   auto result = std::make_unique<TrajectorySeedCollection>();
0189   auto resultTracks = std::make_unique<reco::TrackCollection>();
0190 
0191   const tensorflow::TensorShape input_size_eta({1, 1});
0192   const tensorflow::TensorShape input_size_pt({1, 1});
0193   const tensorflow::TensorShape input_size_cluster({1, jetDimX, jetDimY, Nlayer});
0194   std::vector<std::string> output_names;
0195   output_names.push_back(outputTensorName_[0]);
0196   output_names.push_back(outputTensorName_[1]);
0197 
0198   using namespace edm;
0199   using namespace reco;
0200 
0201   geometry_ = &iSetup.getData(geometryToken_);
0202 
0203   const auto& inputPixelClusters_ = iEvent.get(pixelClusters_);
0204   const auto& vertices = iEvent.get(vertices_);
0205   const auto& cores = iEvent.get(cores_);
0206 
0207   const PixelClusterParameterEstimator* pixelCPE = &iSetup.getData(pixelCPEToken_);
0208 
0209   const TrackerTopology* const tTopo = &iSetup.getData(topoToken_);
0210   auto output = std::make_unique<edmNew::DetSetVector<SiPixelCluster>>();
0211 
0212   for (const auto& jet : cores) {  // jet loop
0213 
0214     if (jet.pt() > ptMin_) {
0215       std::set<unsigned long long int> ids;
0216       const reco::Vertex& jetVertex = vertices[0];
0217 
0218       std::vector<GlobalVector> splitClustDirSet =
0219           splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 2, inputPixelClusters_);
0220       bool l2off = (splitClustDirSet.empty());  // no adc BPIX2
0221       splitClustDirSet = splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 3, inputPixelClusters_);
0222       bool l134off = (splitClustDirSet.empty());  // no adc BPIX1 or BPIX3 or BPIX4
0223       splitClustDirSet = splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 4, inputPixelClusters_);
0224       l134off = (splitClustDirSet.empty() && l134off);
0225       splitClustDirSet = splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 1, inputPixelClusters_);
0226       l134off = (splitClustDirSet.empty() && l134off);
0227 
0228       if (splitClustDirSet.empty()) {  //if layer 1 is broken find direcitons on layer 2
0229         splitClustDirSet = splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 2, inputPixelClusters_);
0230       }
0231       splitClustDirSet.emplace_back(GlobalVector(jet.px(), jet.py(), jet.pz()));
0232       for (int cc = 0; cc < (int)splitClustDirSet.size(); cc++) {
0233         tensorflow::NamedTensorList input_tensors;  //TensorFlow tensors init
0234         input_tensors.resize(3);
0235         input_tensors[0] =
0236             tensorflow::NamedTensor(inputTensorName_[0], tensorflow::Tensor(tensorflow::DT_FLOAT, input_size_eta));
0237         input_tensors[1] =
0238             tensorflow::NamedTensor(inputTensorName_[1], tensorflow::Tensor(tensorflow::DT_FLOAT, input_size_pt));
0239         input_tensors[2] = tensorflow::NamedTensor(inputTensorName_[2],
0240                                                    tensorflow::Tensor(tensorflow::DT_FLOAT, {input_size_cluster}));
0241 
0242         //put all the input tensor to 0
0243         input_tensors[0].second.matrix<float>()(0, 0) = 0.0;
0244         input_tensors[1].second.matrix<float>()(0, 0) = 0.0;
0245         for (int x = 0; x < jetDimX; x++) {
0246           for (int y = 0; y < jetDimY; y++) {
0247             for (int l = 0; l < Nlayer; l++) {
0248               input_tensors[2].second.tensor<float, 4>()(0, x, y, l) = 0.0;
0249             }
0250           }
0251         }  //end of TensorFlow tensors init
0252 
0253         GlobalVector bigClustDir = splitClustDirSet[cc];
0254 
0255         jetEta_ = jet.eta();
0256         jetPt_ = jet.pt();
0257         input_tensors[0].second.matrix<float>()(0, 0) = jet.eta();
0258         input_tensors[1].second.matrix<float>()(0, 0) = jet.pt();
0259 
0260         const GeomDet* globDet = DetectorSelector(
0261             2, jet, bigClustDir, jetVertex, tTopo, inputPixelClusters_);  //det. aligned to bigClustDir on L2
0262 
0263         if (globDet == nullptr)  //no intersection between bigClustDir and pixel detector modules found
0264           continue;
0265 
0266         const GeomDet* goodDet1 = DetectorSelector(1, jet, bigClustDir, jetVertex, tTopo, inputPixelClusters_);
0267         const GeomDet* goodDet3 = DetectorSelector(3, jet, bigClustDir, jetVertex, tTopo, inputPixelClusters_);
0268         const GeomDet* goodDet4 = DetectorSelector(4, jet, bigClustDir, jetVertex, tTopo, inputPixelClusters_);
0269 
0270         for (const auto& detset : inputPixelClusters_) {
0271           const GeomDet* det = geometry_->idToDet(detset.id());
0272 
0273           for (const auto& aCluster : detset) {
0274             det_id_type aClusterID = detset.id();
0275             if (DetId(aClusterID).subdetId() != PixelSubdetector::PixelBarrel)
0276               continue;
0277 
0278             int lay = tTopo->layer(det->geographicalId());
0279 
0280             std::pair<bool, Basic3DVector<float>> interPair =
0281                 findIntersection(bigClustDir, (reco::Candidate::Point)jetVertex.position(), det);
0282             if (interPair.first == false)
0283               continue;
0284             Basic3DVector<float> inter = interPair.second;
0285             auto localInter = det->specificSurface().toLocal((GlobalPoint)inter);
0286 
0287             GlobalPoint pointVertex(jetVertex.position().x(), jetVertex.position().y(), jetVertex.position().z());
0288 
0289             LocalPoint clustPos_local =
0290                 pixelCPE->localParametersV(aCluster, (*geometry_->idToDetUnit(detset.id())))[0].first;
0291 
0292             if (std::abs(clustPos_local.x() - localInter.x()) / pitchX_ <= jetDimX / 2 &&
0293                 std::abs(clustPos_local.y() - localInter.y()) / pitchY_ <=
0294                     jetDimY / 2) {  //used the baricenter, better description maybe useful
0295 
0296               if (det == goodDet1 || det == goodDet3 || det == goodDet4 || det == globDet) {
0297                 fillPixelMatrix(aCluster, lay, localInter, det, input_tensors);
0298               }
0299             }  //cluster in ROI
0300           }  //cluster
0301         }  //detset
0302 
0303         //here the NN produce the seed from the filled input
0304         std::pair<double[jetDimX][jetDimY][Nover][Npar], double[jetDimX][jetDimY][Nover]> seedParamNN =
0305             DeepCoreSeedGenerator::SeedEvaluation(input_tensors, output_names);
0306 
0307         for (int i = 0; i < jetDimX; i++) {
0308           for (int j = 0; j < jetDimY; j++) {
0309             for (int o = 0; o < Nover; o++) {
0310               if (seedParamNN.second[i][j][o] >
0311                   (probThr_ + dth[o] + (l2off ? 0.3 : 0) + (l134off ? dthl[o] : 0))) {  // DeepCore 2.2.1 Threshold
0312                 std::pair<bool, Basic3DVector<float>> interPair =
0313                     findIntersection(bigClustDir, (reco::Candidate::Point)jetVertex.position(), globDet);
0314                 auto localInter = globDet->specificSurface().toLocal((GlobalPoint)interPair.second);
0315 
0316                 int flip = pixelFlipper(globDet);  // 1=not flip, -1=flip
0317                 int nx = i - jetDimX / 2;
0318                 int ny = j - jetDimY / 2;
0319                 nx = flip * nx;
0320                 std::pair<int, int> pixInter = local2Pixel(localInter.x(), localInter.y(), globDet);
0321                 nx = nx + pixInter.first;
0322                 ny = ny + pixInter.second;
0323                 LocalPoint xyLocal = pixel2Local(nx, ny, globDet);
0324 
0325                 double xx = xyLocal.x() + seedParamNN.first[i][j][o][0] * 0.01;  //0.01=internal normalization of NN
0326                 double yy = xyLocal.y() + seedParamNN.first[i][j][o][1] * 0.01;  //0.01=internal normalization of NN
0327                 LocalPoint localSeedPoint = LocalPoint(xx, yy, 0);
0328 
0329                 double track_eta =
0330                     seedParamNN.first[i][j][o][2] * 0.01 + bigClustDir.eta();  //pay attention to this 0.01
0331                 double track_theta = 2 * std::atan(std::exp(-track_eta));
0332                 double track_phi =
0333                     seedParamNN.first[i][j][o][3] * 0.01 + bigClustDir.phi();  //pay attention to this 0.01
0334                 double pt = seedParamNN.first[i][j][o][4];
0335                 double normdirR = pt / sin(track_theta);
0336 
0337                 const GlobalVector globSeedDir(
0338                     GlobalVector::Polar(Geom::Theta<double>(track_theta), Geom::Phi<double>(track_phi), normdirR));
0339                 LocalVector localSeedDir = globDet->surface().toLocal(globSeedDir);
0340                 uint64_t seedid = (uint64_t(xx * 200.) << 0) + (uint64_t(yy * 200.) << 16) +
0341                                   (uint64_t(track_eta * 400.) << 32) + (uint64_t(track_phi * 400.) << 48);
0342                 if (ids.count(seedid) != 0) {
0343                   continue;
0344                 }
0345                 //to disable jetcore iteration comment from here to probability block end--->
0346                 //seeing iteration skipped, useful to total eff and FakeRate comparison
0347                 ids.insert(seedid);
0348 
0349                 //Covariance matrix, the hadrcoded variances = NN residuals width
0350                 //(see https://twiki.cern.ch/twiki/bin/view/CMSPublic/NNJetCoreAtCtD2019)
0351                 float em[15] = {0,
0352                                 0,
0353                                 0,
0354                                 0,
0355                                 0,
0356                                 0,
0357                                 0,
0358                                 0,
0359                                 0,
0360                                 0,
0361                                 0,
0362                                 0,
0363                                 0,
0364                                 0,
0365                                 0};   // (see LocalTrajectoryError for details), order as follow:
0366                 em[0] = 0.15 * 0.15;  // q/pt
0367                 em[2] = 0.5e-5;       // dxdz
0368                 em[5] = 0.5e-5;       // dydz
0369                 em[9] = 2e-5;         // x
0370                 em[14] = 2e-5;        // y
0371                 long int detId = globDet->geographicalId();
0372                 LocalTrajectoryParameters localParam(localSeedPoint, localSeedDir, TrackCharge(1));
0373                 result->emplace_back(TrajectorySeed(PTrajectoryStateOnDet(localParam, pt, em, detId, /*surfaceSide*/ 0),
0374                                                     edm::OwnVector<TrackingRecHit>(),
0375                                                     PropagationDirection::alongMomentum));
0376 
0377                 GlobalPoint globalSeedPoint = globDet->surface().toGlobal(localSeedPoint);
0378                 reco::Track::CovarianceMatrix mm;
0379                 resultTracks->emplace_back(
0380                     reco::Track(1,
0381                                 1,
0382                                 reco::Track::Point(globalSeedPoint.x(), globalSeedPoint.y(), globalSeedPoint.z()),
0383                                 reco::Track::Vector(globSeedDir.x(), globSeedDir.y(), globSeedDir.z()),
0384                                 1,
0385                                 mm));
0386               }  //probability
0387             }
0388           }
0389         }
0390       }  //bigcluster
0391     }  //jet > pt
0392   }  //jet
0393   iEvent.put(std::move(result));
0394   iEvent.put(std::move(resultTracks));
0395 }
0396 
0397 std::pair<bool, Basic3DVector<float>> DeepCoreSeedGenerator::findIntersection(const GlobalVector& dir,
0398                                                                               const reco::Candidate::Point& vertex,
0399                                                                               const GeomDet* det) {
0400   StraightLinePlaneCrossing vertexPlane(Basic3DVector<float>(vertex.x(), vertex.y(), vertex.z()),
0401                                         Basic3DVector<float>(dir.x(), dir.y(), dir.z()));
0402 
0403   std::pair<bool, Basic3DVector<float>> pos = vertexPlane.position(det->specificSurface());
0404 
0405   return pos;
0406 }
0407 
0408 std::pair<int, int> DeepCoreSeedGenerator::local2Pixel(double locX, double locY, const GeomDet* det) {
0409   LocalPoint locXY(locX, locY);
0410   float pixX = (dynamic_cast<const PixelGeomDetUnit*>(det))->specificTopology().pixel(locXY).first;
0411   float pixY = (dynamic_cast<const PixelGeomDetUnit*>(det))->specificTopology().pixel(locXY).second;
0412   std::pair<int, int> out(pixX, pixY);
0413   return out;
0414 }
0415 
0416 LocalPoint DeepCoreSeedGenerator::pixel2Local(int pixX, int pixY, const GeomDet* det) {
0417   float locX = (dynamic_cast<const PixelGeomDetUnit*>(det))->specificTopology().localX(pixX);
0418   float locY = (dynamic_cast<const PixelGeomDetUnit*>(det))->specificTopology().localY(pixY);
0419   LocalPoint locXY(locX, locY);
0420   return locXY;
0421 }
0422 
0423 int DeepCoreSeedGenerator::pixelFlipper(const GeomDet* det) {
0424   int out = 1;
0425   LocalVector locZdir(0, 0, 1);
0426   GlobalVector globZdir = det->specificSurface().toGlobal(locZdir);
0427   const GlobalPoint& globDetCenter = det->position();
0428   float direction =
0429       globZdir.x() * globDetCenter.x() + globZdir.y() * globDetCenter.y() + globZdir.z() * globDetCenter.z();
0430   if (direction < 0)
0431     out = -1;
0432   return out;
0433 }
0434 
0435 void DeepCoreSeedGenerator::fillPixelMatrix(const SiPixelCluster& cluster,
0436                                             int layer,
0437                                             Point3DBase<float, LocalTag> inter,
0438                                             const GeomDet* det,
0439                                             tensorflow::NamedTensorList input_tensors) {
0440   int flip = pixelFlipper(det);  // 1=not flip, -1=flip
0441 
0442   for (int i = 0; i < cluster.size(); i++) {
0443     SiPixelCluster::Pixel pix = cluster.pixel(i);
0444     std::pair<int, int> pixInter = local2Pixel(inter.x(), inter.y(), det);
0445     int nx = pix.x - pixInter.first;
0446     int ny = pix.y - pixInter.second;
0447     nx = flip * nx;
0448 
0449     if (abs(nx) < jetDimX / 2 && abs(ny) < jetDimY / 2) {
0450       nx = nx + jetDimX / 2;
0451       ny = ny + jetDimY / 2;
0452       //14000 = normalization of ACD counts used in the NN
0453       input_tensors[2].second.tensor<float, 4>()(0, nx, ny, layer - 1) += (pix.adc) / (14000.f);
0454     }
0455   }
0456 }
0457 
0458 std::pair<double[DeepCoreSeedGenerator::jetDimX][DeepCoreSeedGenerator::jetDimY][DeepCoreSeedGenerator::Nover]
0459                 [DeepCoreSeedGenerator::Npar],
0460           double[DeepCoreSeedGenerator::jetDimX][DeepCoreSeedGenerator::jetDimY][DeepCoreSeedGenerator::Nover]>
0461 DeepCoreSeedGenerator::SeedEvaluation(tensorflow::NamedTensorList input_tensors,
0462                                       std::vector<std::string> output_names) {
0463   std::vector<tensorflow::Tensor> outputs;
0464   tensorflow::run(session_, input_tensors, output_names, &outputs);
0465   auto matrix_output_par = outputs.at(0).tensor<float, 5>();
0466   auto matrix_output_prob = outputs.at(1).tensor<float, 5>();
0467 
0468   std::pair<double[jetDimX][jetDimY][Nover][Npar], double[jetDimX][jetDimY][Nover]> output_combined;
0469 
0470   for (int x = 0; x < jetDimX; x++) {
0471     for (int y = 0; y < jetDimY; y++) {
0472       for (int trk = 0; trk < Nover; trk++) {
0473         output_combined.second[x][y][trk] =
0474             matrix_output_prob(0, x, y, trk, 0);  //outputs.at(1).matrix<double>()(0,x,y,trk);
0475 
0476         for (int p = 0; p < Npar; p++) {
0477           output_combined.first[x][y][trk][p] =
0478               matrix_output_par(0, x, y, trk, p);  //outputs.at(0).matrix<double>()(0,x,y,trk,p);
0479         }
0480       }
0481     }
0482   }
0483   return output_combined;
0484 }
0485 
0486 const GeomDet* DeepCoreSeedGenerator::DetectorSelector(int llay,
0487                                                        const reco::Candidate& jet,
0488                                                        GlobalVector jetDir,
0489                                                        const reco::Vertex& jetVertex,
0490                                                        const TrackerTopology* const tTopo,
0491                                                        const edmNew::DetSetVector<SiPixelCluster>& clusters) {
0492   double minDist = 0.0;
0493   GeomDet* output = (GeomDet*)nullptr;
0494   for (const auto& detset : clusters) {
0495     auto aClusterID = detset.id();
0496     if (DetId(aClusterID).subdetId() != 1)
0497       continue;
0498     const GeomDet* det = geometry_->idToDet(aClusterID);
0499     int lay = tTopo->layer(det->geographicalId());
0500     if (lay != llay)
0501       continue;
0502     std::pair<bool, Basic3DVector<float>> interPair =
0503         findIntersection(jetDir, (reco::Candidate::Point)jetVertex.position(), det);
0504     if (interPair.first == false)
0505       continue;
0506     Basic3DVector<float> inter = interPair.second;
0507     auto localInter = det->specificSurface().toLocal((GlobalPoint)inter);
0508     if ((minDist == 0.0 || std::abs(localInter.x()) < minDist) && std::abs(localInter.y()) < 3.35) {
0509       minDist = std::abs(localInter.x());
0510       output = (GeomDet*)det;
0511     }
0512   }  //detset
0513   return output;
0514 }
0515 
0516 std::vector<GlobalVector> DeepCoreSeedGenerator::splittedClusterDirections(
0517     const reco::Candidate& jet,
0518     const TrackerTopology* const tTopo,
0519     const PixelClusterParameterEstimator* pixelCPE,
0520     const reco::Vertex& jetVertex,
0521     int layer,
0522     const edmNew::DetSetVector<SiPixelCluster>& clusters) {
0523   std::vector<GlobalVector> clustDirs;
0524   for (const auto& detset_int : clusters) {
0525     const GeomDet* det_int = geometry_->idToDet(detset_int.id());
0526     int lay = tTopo->layer(det_int->geographicalId());
0527     if (lay != layer)
0528       continue;  //NB: saved bigClusters on all the layers!!
0529     auto detUnit = geometry_->idToDetUnit(detset_int.id());
0530     for (const auto& aCluster : detset_int) {
0531       GlobalPoint clustPos = det_int->surface().toGlobal(pixelCPE->localParametersV(aCluster, (*detUnit))[0].first);
0532       GlobalPoint vertexPos(jetVertex.position().x(), jetVertex.position().y(), jetVertex.position().z());
0533       GlobalVector clusterDir = clustPos - vertexPos;
0534       GlobalVector jetDir(jet.px(), jet.py(), jet.pz());
0535       if (Geom::deltaR(jetDir, clusterDir) < deltaR_) {
0536         clustDirs.emplace_back(clusterDir);
0537       }
0538     }
0539   }
0540   return clustDirs;
0541 }
0542 
0543 std::unique_ptr<tensorflow::SessionCache> DeepCoreSeedGenerator::initializeGlobalCache(
0544     const edm::ParameterSet& iConfig) {
0545   // this method is supposed to create, initialize and return a Tensorflow:SessionCache instance
0546   std::string graphPath = iConfig.getParameter<edm::FileInPath>("weightFile").fullPath();
0547   std::unique_ptr<tensorflow::SessionCache> cache = std::make_unique<tensorflow::SessionCache>(graphPath);
0548   return cache;
0549 }
0550 
0551 void DeepCoreSeedGenerator::globalEndJob(tensorflow::SessionCache* cache) {}
0552 
0553 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0554 void DeepCoreSeedGenerator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0555   edm::ParameterSetDescription desc;
0556   desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
0557   desc.add<edm::InputTag>("pixelClusters", edm::InputTag("siPixelClustersPreSplitting"));
0558   desc.add<edm::InputTag>("cores", edm::InputTag("jetsForCoreTracking"));
0559   desc.add<double>("ptMin", 100);    // the current training uses 500 GeV
0560   desc.add<double>("deltaR", 0.25);  // the current training uses 0.1
0561   desc.add<double>("chargeFractionMin", 18000.0);
0562   desc.add<double>("centralMIPCharge", 2);
0563   desc.add<std::string>("pixelCPE", "PixelCPEGeneric");
0564   desc.add<edm::FileInPath>(
0565       "weightFile",
0566       edm::FileInPath("RecoTracker/TkSeedGenerator/data/DeepCore/DeepCoreSeedGenerator_TrainedModel_barrel_2023.pb"));
0567   desc.add<std::vector<std::string>>("inputTensorName", {"input_1", "input_2", "input_3"});
0568   desc.add<std::vector<std::string>>("outputTensorName", {"output_node0", "output_node1"});
0569   desc.add<double>("probThr", 0.7);  // DeepCore 2.2.1 baseline threshold
0570   descriptions.add("deepCoreSeedGenerator", desc);
0571 }
0572 
0573 DEFINE_FWK_MODULE(DeepCoreSeedGenerator);
[ 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!