From edd48ca5a6d71bc4a668bcba01d5ae9a2cb0332d Mon Sep 17 00:00:00 2001
From: Hichem Bouchamaoui <bouchamaouihichem@gmail.com>
Date: Mon, 18 Sep 2023 17:06:47 -0500
Subject: [PATCH] Updating DeepCore Seed Generatorto to use DeepCore model
 2.2.1 (2023)

---
 .../plugins/DeepCoreSeedGenerator.cc          | 40 +++++++++++++++----
 1 file changed, 33 insertions(+), 7 deletions(-)

diff --git a/RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.cc b/RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.cc
index e8d4c68fc3ae4..66b5d8e656cac 100644
--- a/RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.cc
+++ b/RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.cc
@@ -95,6 +95,10 @@ class DeepCoreSeedGenerator : public edm::stream::EDProducer<edm::GlobalCache<te
   static constexpr int Nover = 3;     //Max number of tracks recorded per pixel
   static constexpr int Npar = 5;      //Number of track parameter
 
+  // DeepCore 2.2.1 thresholds delta
+  double dth[3] = {0.0, 0.15, 0.3};
+  double dthl[3] = {0.1, 0.05, 0};
+
 private:
   void produce(edm::Event&, const edm::EventSetup&) override;
 
@@ -210,10 +214,25 @@ void DeepCoreSeedGenerator::produce(edm::Event& iEvent, const edm::EventSetup& i
     if (jet.pt() > ptMin_) {
       std::set<unsigned long long int> ids;
       const reco::Vertex& jetVertex = vertices[0];
-
+      // DeepCore 1.0:
+      /*
       std::vector<GlobalVector> splitClustDirSet =
-          splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 1, inputPixelClusters_);
+                  splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 1, inputPixelClusters_);
       bool l2off = (splitClustDirSet.empty());
+
+      //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+      // DeepCore 2.1.3 Threshold checks:
+      */
+      std::vector<GlobalVector> splitClustDirSet =
+          splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 2, inputPixelClusters_);
+      bool l2off = (splitClustDirSet.empty());  // no adc BPIX2
+      splitClustDirSet = splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 3, inputPixelClusters_);
+      bool l134off = (splitClustDirSet.empty());  // no adc BPIX1 or BPIX3 or BPIX4
+      splitClustDirSet = splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 4, inputPixelClusters_);
+      l134off = (splitClustDirSet.empty() && l134off);
+      splitClustDirSet = splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 1, inputPixelClusters_);
+      l134off = (splitClustDirSet.empty() && l134off);
+
       if (splitClustDirSet.empty()) {  //if layer 1 is broken find direcitons on layer 2
         splitClustDirSet = splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 2, inputPixelClusters_);
       }
@@ -296,7 +315,9 @@ void DeepCoreSeedGenerator::produce(edm::Event& iEvent, const edm::EventSetup& i
         for (int i = 0; i < jetDimX; i++) {
           for (int j = 0; j < jetDimY; j++) {
             for (int o = 0; o < Nover; o++) {
-              if (seedParamNN.second[i][j][o] > (probThr_ - o * 0.1 - (l2off ? 0.35 : 0))) {
+              //  if (seedParamNN.second[i][j][o] > (probThr_ - o * 0.1 - (l2off ? 0.35 : 0))) {  // DeepCore 1.0 Threshold
+              if (seedParamNN.second[i][j][o] >
+                  (probThr_ + dth[o] + (l2off ? 0.3 : 0) + (l134off ? dthl[o] : 0))) {  // DeepCore 2.2.1 Threshold
                 std::pair<bool, Basic3DVector<float>> interPair =
                     findIntersection(bigClustDir, (reco::Candidate::Point)jetVertex.position(), globDet);
                 auto localInter = globDet->specificSurface().toLocal((GlobalPoint)interPair.second);
@@ -319,8 +340,9 @@ void DeepCoreSeedGenerator::produce(edm::Event& iEvent, const edm::EventSetup& i
                 double track_theta = 2 * std::atan(std::exp(-track_eta));
                 double track_phi =
                     seedParamNN.first[i][j][o][3] * 0.01 + bigClustDir.phi();  //pay attention to this 0.01
-
-                double pt = 1. / seedParamNN.first[i][j][o][4];
+                // DeepCore 1.0 predicts 1/pt instead of pt, while DeepCore 2.2.1 directly predicts pt
+                // double pt = 1. / seedParamNN.first[i][j][o][4];
+                double pt = seedParamNN.first[i][j][o][4];
                 double normdirR = pt / sin(track_theta);
 
                 const GlobalVector globSeedDir(
@@ -552,10 +574,14 @@ void DeepCoreSeedGenerator::fillDescriptions(edm::ConfigurationDescriptions& des
   desc.add<std::string>("pixelCPE", "PixelCPEGeneric");
   desc.add<edm::FileInPath>(
       "weightFile",
-      edm::FileInPath("RecoTracker/TkSeedGenerator/data/DeepCore/DeepCoreSeedGenerator_TrainedModel_barrel_2017.pb"));
+      // DeepCore 1.0
+      // edm::FileInPath("RecoTracker/TkSeedGenerator/data/DeepCore/DeepCoreSeedGenerator_TrainedModel_barrel_2017.pb"));
+      // DeepCore 2.2.1
+      edm::FileInPath("RecoTracker/TkSeedGenerator/data/DeepCore/DeepCoreSeedGenerator_TrainedModel_barrel_2023.pb"));
   desc.add<std::vector<std::string>>("inputTensorName", {"input_1", "input_2", "input_3"});
   desc.add<std::vector<std::string>>("outputTensorName", {"output_node0", "output_node1"});
-  desc.add<double>("probThr", 0.85);
+  //  desc.add<double>("probThr", 0.85); // DeepCore 1.0
+  desc.add<double>("probThr", 0.7);  // DeepCore 2.2.1 Threshold
   descriptions.add("deepCoreSeedGenerator", desc);
 }