Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-18 02:30:17

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