TkSeedGenerator/plugins/DeepCoreSeedGenerator.cc

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