Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-29 07:47:35

0001 /*
0002  * \class L2TauTagProducer
0003  *
0004  * L2Tau identification using Convolutional NN.
0005  *
0006  * \author Valeria D'Amante, Università di Siena and INFN Pisa
0007  *         Konstantin Androsov, EPFL and ETHZ
0008 */
0009 #include <memory>
0010 #include <boost/property_tree/json_parser.hpp>
0011 #include <cmath>
0012 #include "FWCore/Framework/interface/stream/EDProducer.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "DataFormats/Math/interface/deltaR.h"
0019 #include "DataFormats/Common/interface/Handle.h"
0020 #include "FWCore/Utilities/interface/InputTag.h"
0021 #include "FWCore/Utilities/interface/isFinite.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0024 #include "PhysicsTools/TensorFlow/interface/TensorFlow.h"
0025 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0026 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0027 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0028 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0029 #include "DataFormats/CaloRecHit/interface/CaloRecHit.h"
0030 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0031 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0032 #include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
0033 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0034 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
0035 #include "DataFormats/HcalRecHit/interface/HcalRecHitDefs.h"
0036 #include "DataFormats/HcalRecHit/interface/HFRecHit.h"
0037 #include "DataFormats/HcalRecHit/interface/HORecHit.h"
0038 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
0039 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0040 #include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
0041 #include "RecoPixelVertexing/PixelTrackFitting/interface/FitUtils.h"
0042 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0043 #include "DataFormats/TrackReco/interface/HitPattern.h"
0044 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h"
0045 #include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
0046 #include "DataFormats/GeometrySurface/interface/Plane.h"
0047 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0048 #include "CUDADataFormats/Track/interface/PixelTrackHeterogeneous.h"
0049 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0050 #include "CUDADataFormats/SiPixelCluster/interface/gpuClusteringConstants.h"
0051 #include "CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h"
0052 #include "CUDADataFormats/Vertex/interface/ZVertexSoA.h"
0053 #include "CUDADataFormats/Vertex/interface/ZVertexHeterogeneous.h"
0054 
0055 namespace L2TauTagNNv1 {
0056   constexpr int nCellEta = 5;
0057   constexpr int nCellPhi = 5;
0058   constexpr int nVars = 31;
0059   constexpr float dR_max = 0.5;
0060   enum class NNInputs {
0061     nVertices = 0,
0062     l1Tau_pt,
0063     l1Tau_eta,
0064     l1Tau_hwIso,
0065     EcalEnergySum,
0066     EcalSize,
0067     EcalEnergyStdDev,
0068     EcalDeltaEta,
0069     EcalDeltaPhi,
0070     EcalChi2,
0071     EcalEnergySumForPositiveChi2,
0072     EcalSizeForPositiveChi2,
0073     HcalEnergySum,
0074     HcalSize,
0075     HcalEnergyStdDev,
0076     HcalDeltaEta,
0077     HcalDeltaPhi,
0078     HcalChi2,
0079     HcalEnergySumForPositiveChi2,
0080     HcalSizeForPositiveChi2,
0081     PatatrackPtSum,
0082     PatatrackSize,
0083     PatatrackSizeWithVertex,
0084     PatatrackPtSumWithVertex,
0085     PatatrackChargeSum,
0086     PatatrackDeltaEta,
0087     PatatrackDeltaPhi,
0088     PatatrackChi2OverNdof,
0089     PatatrackNdof,
0090     PatatrackDxy,
0091     PatatrackDz
0092   };
0093 
0094   const std::map<NNInputs, std::string> varNameMap = {
0095       {NNInputs::nVertices, "nVertices"},
0096       {NNInputs::l1Tau_pt, "l1Tau_pt"},
0097       {NNInputs::l1Tau_eta, "l1Tau_eta"},
0098       {NNInputs::l1Tau_hwIso, "l1Tau_hwIso"},
0099       {NNInputs::EcalEnergySum, "EcalEnergySum"},
0100       {NNInputs::EcalSize, "EcalSize"},
0101       {NNInputs::EcalEnergyStdDev, "EcalEnergyStdDev"},
0102       {NNInputs::EcalDeltaEta, "EcalDeltaEta"},
0103       {NNInputs::EcalDeltaPhi, "EcalDeltaPhi"},
0104       {NNInputs::EcalChi2, "EcalChi2"},
0105       {NNInputs::EcalEnergySumForPositiveChi2, "EcalEnergySumForPositiveChi2"},
0106       {NNInputs::EcalSizeForPositiveChi2, "EcalSizeForPositiveChi2"},
0107       {NNInputs::HcalEnergySum, "HcalEnergySum"},
0108       {NNInputs::HcalSize, "HcalSize"},
0109       {NNInputs::HcalEnergyStdDev, "HcalEnergyStdDev"},
0110       {NNInputs::HcalDeltaEta, "HcalDeltaEta"},
0111       {NNInputs::HcalDeltaPhi, "HcalDeltaPhi"},
0112       {NNInputs::HcalChi2, "HcalChi2"},
0113       {NNInputs::HcalEnergySumForPositiveChi2, "HcalEnergySumForPositiveChi2"},
0114       {NNInputs::HcalSizeForPositiveChi2, "HcalSizeForPositiveChi2"},
0115       {NNInputs::PatatrackPtSum, "PatatrackPtSum"},
0116       {NNInputs::PatatrackSize, "PatatrackSize"},
0117       {NNInputs::PatatrackSizeWithVertex, "PatatrackSizeWithVertex"},
0118       {NNInputs::PatatrackPtSumWithVertex, "PatatrackPtSumWithVertex"},
0119       {NNInputs::PatatrackChargeSum, "PatatrackChargeSum"},
0120       {NNInputs::PatatrackDeltaEta, "PatatrackDeltaEta"},
0121       {NNInputs::PatatrackDeltaPhi, "PatatrackDeltaPhi"},
0122       {NNInputs::PatatrackChi2OverNdof, "PatatrackChi2OverNdof"},
0123       {NNInputs::PatatrackNdof, "PatatrackNdof"},
0124       {NNInputs::PatatrackDxy, "PatatrackDxy"},
0125       {NNInputs::PatatrackDz, "PatatrackDz"}};
0126 }  // namespace L2TauTagNNv1
0127 namespace {
0128   inline float& getCellImpl(
0129       tensorflow::Tensor& cellGridMatrix, int tau_idx, int phi_idx, int eta_idx, L2TauTagNNv1::NNInputs NNInput_idx) {
0130     return cellGridMatrix.tensor<float, 4>()(tau_idx, phi_idx, eta_idx, static_cast<int>(NNInput_idx));
0131   }
0132 }  // namespace
0133 struct normDictElement {
0134   float mean;
0135   float std;
0136   float min;
0137   float max;
0138 };
0139 
0140 struct L2TauNNProducerCacheData {
0141   L2TauNNProducerCacheData() : graphDef(nullptr), session(nullptr) {}
0142   tensorflow::GraphDef* graphDef;
0143   tensorflow::Session* session;
0144   std::vector<normDictElement> normVec;
0145 };
0146 
0147 class L2TauNNProducer : public edm::stream::EDProducer<edm::GlobalCache<L2TauNNProducerCacheData>> {
0148 public:
0149   struct caloRecHitCollections {
0150     const HBHERecHitCollection* hbhe;
0151     const HORecHitCollection* ho;
0152     const EcalRecHitCollection* eb;
0153     const EcalRecHitCollection* ee;
0154     const CaloGeometry* geometry;
0155   };
0156 
0157   struct InputDescTau {
0158     std::string CollectionName;
0159     edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> inputToken_;
0160   };
0161 
0162   static constexpr float dR2_max = L2TauTagNNv1::dR_max * L2TauTagNNv1::dR_max;
0163   static constexpr float dEta_width = 2 * L2TauTagNNv1::dR_max / static_cast<float>(L2TauTagNNv1::nCellEta);
0164   static constexpr float dPhi_width = 2 * L2TauTagNNv1::dR_max / static_cast<float>(L2TauTagNNv1::nCellPhi);
0165 
0166   explicit L2TauNNProducer(const edm::ParameterSet&, const L2TauNNProducerCacheData*);
0167   static void fillDescriptions(edm::ConfigurationDescriptions&);
0168   static std::unique_ptr<L2TauNNProducerCacheData> initializeGlobalCache(const edm::ParameterSet&);
0169   static void globalEndJob(L2TauNNProducerCacheData*);
0170 
0171 private:
0172   void checknan(tensorflow::Tensor& tensor, int debugLevel);
0173   void standardizeTensor(tensorflow::Tensor& tensor);
0174   std::vector<float> getTauScore(const tensorflow::Tensor& cellGridMatrix);
0175   void produce(edm::Event& event, const edm::EventSetup& eventsetup) override;
0176   void fillL1TauVars(tensorflow::Tensor& cellGridMatrix, const std::vector<l1t::TauRef>& allTaus);
0177   void fillCaloRecHits(tensorflow::Tensor& cellGridMatrix,
0178                        const std::vector<l1t::TauRef>& allTaus,
0179                        const caloRecHitCollections& caloRecHits);
0180   void fillPatatracks(tensorflow::Tensor& cellGridMatrix,
0181                       const std::vector<l1t::TauRef>& allTaus,
0182                       const pixelTrack::TrackSoA& patatracks_tsoa,
0183                       const ZVertexSoA& patavtx_soa,
0184                       const reco::BeamSpot& beamspot,
0185                       const MagneticField* magfi);
0186   void selectGoodTracksAndVertices(const ZVertexSoA& patavtx_soa,
0187                                    const pixelTrack::TrackSoA& patatracks_tsoa,
0188                                    std::vector<int>& trkGood,
0189                                    std::vector<int>& vtxGood);
0190   std::pair<float, float> impactParameter(int it,
0191                                           const pixelTrack::TrackSoA& patatracks_tsoa,
0192                                           float patatrackPhi,
0193                                           const reco::BeamSpot& beamspot,
0194                                           const MagneticField* magfi);
0195   template <typename VPos, typename LVec>
0196   std::tuple<float, float, int, int> getEtaPhiIndices(const VPos& position, const LVec& tau_p4);
0197   template <typename LVec>
0198   std::tuple<float, float, int, int> getEtaPhiIndices(float eta, float phi, const LVec& tau_p4);
0199 
0200 private:
0201   const int debugLevel_;
0202   const edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> tauTriggerToken_;
0203   std::vector<InputDescTau> L1TauDesc_;
0204   const edm::EDGetTokenT<HBHERecHitCollection> hbheToken_;
0205   const edm::EDGetTokenT<HORecHitCollection> hoToken_;
0206   const edm::EDGetTokenT<EcalRecHitCollection> ebToken_;
0207   const edm::EDGetTokenT<EcalRecHitCollection> eeToken_;
0208   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometryToken_;
0209   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> bFieldToken_;
0210   const edm::EDGetTokenT<ZVertexHeterogeneous> pataVerticesToken_;
0211   const edm::EDGetTokenT<PixelTrackHeterogeneous> pataTracksToken_;
0212   const edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
0213   const unsigned int maxVtx_;
0214   const float fractionSumPt2_;
0215   const float minSumPt2_;
0216   const float trackPtMin_;
0217   const float trackPtMax_;
0218   const float trackChi2Max_;
0219   std::string inputTensorName_;
0220   std::string outputTensorName_;
0221   const L2TauNNProducerCacheData* L2cacheData_;
0222 };
0223 
0224 std::unique_ptr<L2TauNNProducerCacheData> L2TauNNProducer::initializeGlobalCache(const edm::ParameterSet& cfg) {
0225   std::unique_ptr<L2TauNNProducerCacheData> cacheData = std::make_unique<L2TauNNProducerCacheData>();
0226   cacheData->normVec.reserve(L2TauTagNNv1::nVars);
0227 
0228   auto const graphPath = edm::FileInPath(cfg.getParameter<std::string>("graphPath")).fullPath();
0229 
0230   cacheData->graphDef = tensorflow::loadGraphDef(graphPath);
0231   cacheData->session = tensorflow::createSession(cacheData->graphDef);
0232 
0233   tensorflow::setLogging("2");
0234 
0235   boost::property_tree::ptree loadPtreeRoot;
0236   auto const normalizationDict = edm::FileInPath(cfg.getParameter<std::string>("normalizationDict")).fullPath();
0237   boost::property_tree::read_json(normalizationDict, loadPtreeRoot);
0238   for (const auto& [key, val] : L2TauTagNNv1::varNameMap) {
0239     boost::property_tree::ptree var = loadPtreeRoot.get_child(val);
0240     normDictElement current_element;
0241     current_element.mean = var.get_child("mean").get_value<float>();
0242     current_element.std = var.get_child("std").get_value<float>();
0243     current_element.min = var.get_child("min").get_value<float>();
0244     current_element.max = var.get_child("max").get_value<float>();
0245     cacheData->normVec.push_back(current_element);
0246   }
0247   return cacheData;
0248 }
0249 void L2TauNNProducer::globalEndJob(L2TauNNProducerCacheData* cacheData) {
0250   if (cacheData->graphDef != nullptr) {
0251     delete cacheData->graphDef;
0252   }
0253   tensorflow::closeSession(cacheData->session);
0254 }
0255 void L2TauNNProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0256   edm::ParameterSetDescription desc;
0257   desc.add<int>("debugLevel", 0)->setComment("set debug level for printing out info");
0258   edm::ParameterSetDescription l1TausPset;
0259   l1TausPset.add<std::string>("L1CollectionName", "DoubleTau")->setComment("Name of collections");
0260   l1TausPset.add<edm::InputTag>("L1TauTrigger", edm::InputTag("hltL1sDoubleTauBigOR"))
0261       ->setComment("Which trigger should the L1 Taus collection pass");
0262   edm::ParameterSet l1TausPSetDefault;
0263   l1TausPSetDefault.addParameter<std::string>("L1CollectionName", "DoubleTau");
0264   l1TausPSetDefault.addParameter<edm::InputTag>("L1TauTrigger", edm::InputTag("hltL1sDoubleTauBigOR"));
0265   desc.addVPSet("L1Taus", l1TausPset, {l1TausPSetDefault});
0266   desc.add<edm::InputTag>("hbheInput", edm::InputTag("hltHbhereco"))->setComment("HBHE recHit collection");
0267   desc.add<edm::InputTag>("hoInput", edm::InputTag("hltHoreco"))->setComment("HO recHit Collection");
0268   desc.add<edm::InputTag>("ebInput", edm::InputTag("hltEcalRecHit:EcalRecHitsEB"))->setComment("EB recHit Collection");
0269   desc.add<edm::InputTag>("eeInput", edm::InputTag("hltEcalRecHit:EcalRecHitsEE"))->setComment("EE recHit Collection");
0270   desc.add<edm::InputTag>("pataVertices", edm::InputTag("hltPixelVerticesSoA"))
0271       ->setComment("patatrack vertices collection");
0272   desc.add<edm::InputTag>("pataTracks", edm::InputTag("hltPixelTracksSoA"))->setComment("patatrack collection");
0273   desc.add<edm::InputTag>("BeamSpot", edm::InputTag("hltOnlineBeamSpot"))->setComment("BeamSpot Collection");
0274   desc.add<uint>("maxVtx", 100)->setComment("max output collection size (number of accepted vertices)");
0275   desc.add<double>("fractionSumPt2", 0.3)->setComment("threshold on sumPt2 fraction of the leading vertex");
0276   desc.add<double>("minSumPt2", 0.)->setComment("min sumPt2");
0277   desc.add<double>("track_pt_min", 1.0)->setComment("min track p_T");
0278   desc.add<double>("track_pt_max", 10.0)->setComment("max track p_T");
0279   desc.add<double>("track_chi2_max", 99999.)->setComment("max track chi2");
0280   desc.add<std::string>("graphPath", "RecoTauTag/TrainingFiles/data/L2TauNNTag/L2TauTag_Run3v1.pb")
0281       ->setComment("path to the saved CNN");
0282   desc.add<std::string>("normalizationDict", "RecoTauTag/TrainingFiles/data/L2TauNNTag/NormalizationDict.json")
0283       ->setComment("path to the dictionary for variable standardization");
0284   descriptions.addWithDefaultLabel(desc);
0285 }
0286 
0287 L2TauNNProducer::L2TauNNProducer(const edm::ParameterSet& cfg, const L2TauNNProducerCacheData* cacheData)
0288     : debugLevel_(cfg.getParameter<int>("debugLevel")),
0289       hbheToken_(consumes<HBHERecHitCollection>(cfg.getParameter<edm::InputTag>("hbheInput"))),
0290       hoToken_(consumes<HORecHitCollection>(cfg.getParameter<edm::InputTag>("hoInput"))),
0291       ebToken_(consumes<EcalRecHitCollection>(cfg.getParameter<edm::InputTag>("ebInput"))),
0292       eeToken_(consumes<EcalRecHitCollection>(cfg.getParameter<edm::InputTag>("eeInput"))),
0293       geometryToken_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0294       bFieldToken_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0295       pataVerticesToken_(consumes<ZVertexHeterogeneous>(cfg.getParameter<edm::InputTag>("pataVertices"))),
0296       pataTracksToken_(consumes<PixelTrackHeterogeneous>(cfg.getParameter<edm::InputTag>("pataTracks"))),
0297       beamSpotToken_(consumes<reco::BeamSpot>(cfg.getParameter<edm::InputTag>("BeamSpot"))),
0298       maxVtx_(cfg.getParameter<uint>("maxVtx")),
0299       fractionSumPt2_(cfg.getParameter<double>("fractionSumPt2")),
0300       minSumPt2_(cfg.getParameter<double>("minSumPt2")),
0301       trackPtMin_(cfg.getParameter<double>("track_pt_min")),
0302       trackPtMax_(cfg.getParameter<double>("track_pt_max")),
0303       trackChi2Max_(cfg.getParameter<double>("track_chi2_max")) {
0304   if (cacheData->graphDef == nullptr) {
0305     throw cms::Exception("InvalidCacheData") << "Invalid Cache Data.";
0306   }
0307   inputTensorName_ = cacheData->graphDef->node(0).name();
0308   outputTensorName_ = cacheData->graphDef->node(cacheData->graphDef->node_size() - 1).name();
0309   L2cacheData_ = cacheData;
0310   std::vector<edm::ParameterSet> L1TauCollections = cfg.getParameter<std::vector<edm::ParameterSet>>("L1Taus");
0311   L1TauDesc_.reserve(L1TauCollections.size());
0312   for (const auto& l1TauInput : L1TauCollections) {
0313     InputDescTau toInsert;
0314     toInsert.CollectionName = l1TauInput.getParameter<std::string>("L1CollectionName");
0315     toInsert.inputToken_ =
0316         consumes<trigger::TriggerFilterObjectWithRefs>(l1TauInput.getParameter<edm::InputTag>("L1TauTrigger"));
0317     L1TauDesc_.push_back(toInsert);
0318   }
0319   for (const auto& desc : L1TauDesc_)
0320     produces<std::vector<float>>(desc.CollectionName);
0321 }
0322 
0323 void L2TauNNProducer::checknan(tensorflow::Tensor& tensor, int debugLevel) {
0324   using NNInputs = L2TauTagNNv1::NNInputs;
0325   std::vector<int> tensor_shape(tensor.shape().dims());
0326   for (int d = 0; d < tensor.shape().dims(); d++) {
0327     tensor_shape.at(d) = tensor.shape().dim_size(d);
0328   }
0329   if (tensor_shape.size() != 4) {
0330     throw cms::Exception("InvalidTensor") << "Tensor shape does not have 4 dimensions!";
0331   }
0332   for (int tau_idx = 0; tau_idx < tensor_shape.at(0); tau_idx++) {
0333     for (int phi_idx = 0; phi_idx < tensor_shape.at(1); phi_idx++) {
0334       for (int eta_idx = 0; eta_idx < tensor_shape.at(2); eta_idx++) {
0335         for (int var_idx = 0; var_idx < tensor_shape.at(3); var_idx++) {
0336           auto getCell = [&](NNInputs input) -> float& {
0337             return getCellImpl(tensor, tau_idx, phi_idx, eta_idx, input);
0338           };
0339           auto nonstd_var = getCell(static_cast<NNInputs>(var_idx));
0340           if (edm::isNotFinite(nonstd_var)) {
0341             edm::LogWarning("InputVar") << "var is nan \nvar name= "
0342                                         << L2TauTagNNv1::varNameMap.at(static_cast<L2TauTagNNv1::NNInputs>(var_idx))
0343                                         << "\t var_idx = " << var_idx << "\t eta_idx = " << eta_idx
0344                                         << "\t phi_idx = " << phi_idx << "\t tau_idx = " << tau_idx;
0345             if (debugLevel > 2) {
0346               edm::LogWarning("InputVar") << "other vars in same cell \n";
0347               if (var_idx + 1 < tensor_shape.at(3))
0348                 edm::LogWarning("InputVar") << L2TauTagNNv1::varNameMap.at(static_cast<NNInputs>(var_idx + 1))
0349                                             << "\t = " << getCell(static_cast<NNInputs>(var_idx + 1));
0350               if (var_idx + 2 < tensor_shape.at(3))
0351                 edm::LogWarning("InputVar") << L2TauTagNNv1::varNameMap.at(static_cast<NNInputs>(var_idx + 2))
0352                                             << "\t = " << getCell(static_cast<NNInputs>(var_idx + 2));
0353               if (var_idx + 3 < tensor_shape.at(3))
0354                 edm::LogWarning("InputVar") << L2TauTagNNv1::varNameMap.at(static_cast<NNInputs>(var_idx + 3))
0355                                             << "\t = " << getCell(static_cast<NNInputs>(var_idx + 3));
0356               if (var_idx + 4 < tensor_shape.at(3))
0357                 edm::LogWarning("InputVar") << L2TauTagNNv1::varNameMap.at(static_cast<NNInputs>(var_idx + 4))
0358                                             << "\t = " << getCell(static_cast<NNInputs>(var_idx + 4));
0359             }
0360           }
0361         }
0362       }
0363     }
0364   }
0365 }
0366 
0367 void L2TauNNProducer::standardizeTensor(tensorflow::Tensor& tensor) {
0368   using NNInputs = L2TauTagNNv1::NNInputs;
0369   std::vector<int> tensor_shape(tensor.shape().dims());
0370   for (int d = 0; d < tensor.shape().dims(); d++) {
0371     tensor_shape.at(d) = tensor.shape().dim_size(d);
0372   }
0373   if (tensor_shape.size() != 4) {
0374     throw cms::Exception("InvalidTensor") << "Tensor shape does not have 4 dimensions!";
0375   }
0376   for (int tau_idx = 0; tau_idx < tensor_shape.at(0); tau_idx++) {
0377     for (int phi_idx = 0; phi_idx < tensor_shape.at(1); phi_idx++) {
0378       for (int eta_idx = 0; eta_idx < tensor_shape.at(2); eta_idx++) {
0379         for (int var_idx = 0; var_idx < tensor_shape.at(3); var_idx++) {
0380           auto getCell = [&](NNInputs input) -> float& {
0381             return getCellImpl(tensor, tau_idx, phi_idx, eta_idx, input);
0382           };
0383           float mean = L2cacheData_->normVec.at(var_idx).mean;
0384           float std = L2cacheData_->normVec.at(var_idx).std;
0385           float min = L2cacheData_->normVec.at(var_idx).min;
0386           float max = L2cacheData_->normVec.at(var_idx).max;
0387           float nonstd_var = getCell(static_cast<NNInputs>(var_idx));
0388           float std_var = static_cast<float>((nonstd_var - mean) / std);
0389           if (std_var > max) {
0390             std_var = static_cast<float>(max);
0391           } else if (std_var < min) {
0392             std_var = static_cast<float>(min);
0393           }
0394           getCell(static_cast<NNInputs>(var_idx)) = std_var;
0395         }
0396       }
0397     }
0398   }
0399 }
0400 
0401 void L2TauNNProducer::fillL1TauVars(tensorflow::Tensor& cellGridMatrix, const std::vector<l1t::TauRef>& allTaus) {
0402   using NNInputs = L2TauTagNNv1::NNInputs;
0403 
0404   const int nTaus = allTaus.size();
0405   for (int tau_idx = 0; tau_idx < nTaus; tau_idx++) {
0406     for (int eta_idx = 0; eta_idx < L2TauTagNNv1::nCellEta; eta_idx++) {
0407       for (int phi_idx = 0; phi_idx < L2TauTagNNv1::nCellPhi; phi_idx++) {
0408         auto getCell = [&](NNInputs input) -> float& {
0409           return getCellImpl(cellGridMatrix, tau_idx, phi_idx, eta_idx, input);
0410         };
0411         getCell(NNInputs::l1Tau_pt) = allTaus[tau_idx]->pt();
0412         getCell(NNInputs::l1Tau_eta) = allTaus[tau_idx]->eta();
0413         getCell(NNInputs::l1Tau_hwIso) = allTaus[tau_idx]->hwIso();
0414       }
0415     }
0416   }
0417 }
0418 
0419 template <typename LVec>
0420 std::tuple<float, float, int, int> L2TauNNProducer::getEtaPhiIndices(float eta, float phi, const LVec& tau_p4) {
0421   const float deta = eta - tau_p4.eta();
0422   const float dphi = reco::deltaPhi(phi, tau_p4.phi());
0423   const int eta_idx = static_cast<int>(floor((deta + L2TauTagNNv1::dR_max) / dEta_width));
0424   const int phi_idx = static_cast<int>(floor((dphi + L2TauTagNNv1::dR_max) / dPhi_width));
0425   return std::make_tuple(deta, dphi, eta_idx, phi_idx);
0426 }
0427 
0428 template <typename VPos, typename LVec>
0429 std::tuple<float, float, int, int> L2TauNNProducer::getEtaPhiIndices(const VPos& position, const LVec& tau_p4) {
0430   return getEtaPhiIndices(position.eta(), position.phi(), tau_p4);
0431 }
0432 
0433 void L2TauNNProducer::fillCaloRecHits(tensorflow::Tensor& cellGridMatrix,
0434                                       const std::vector<l1t::TauRef>& allTaus,
0435                                       const caloRecHitCollections& caloRecHits) {
0436   using NNInputs = L2TauTagNNv1::NNInputs;
0437 
0438   const int nTaus = allTaus.size();
0439   float deta, dphi;
0440   int eta_idx = 0;
0441   int phi_idx = 0;
0442   int tau_idx = 0;
0443 
0444   auto getCell = [&](NNInputs input) -> float& {
0445     return getCellImpl(cellGridMatrix, tau_idx, phi_idx, eta_idx, input);
0446   };
0447   for (tau_idx = 0; tau_idx < nTaus; tau_idx++) {
0448     // calorechit_EE
0449     for (const auto& caloRecHit_ee : *caloRecHits.ee) {
0450       if (caloRecHit_ee.energy() <= 0)
0451         continue;
0452       const auto& position = caloRecHits.geometry->getGeometry(caloRecHit_ee.id())->getPosition();
0453       const float eeCalEn = caloRecHit_ee.energy();
0454       const float eeCalChi2 = caloRecHit_ee.chi2();
0455       if (reco::deltaR2(position, allTaus[tau_idx]->polarP4()) < dR2_max) {
0456         std::tie(deta, dphi, eta_idx, phi_idx) = getEtaPhiIndices(position, allTaus[tau_idx]->polarP4());
0457         getCell(NNInputs::EcalEnergySum) += eeCalEn;
0458         getCell(NNInputs::EcalSize) += 1.;
0459         getCell(NNInputs::EcalEnergyStdDev) += eeCalEn * eeCalEn;
0460         getCell(NNInputs::EcalDeltaEta) += deta * eeCalEn;
0461         getCell(NNInputs::EcalDeltaPhi) += dphi * eeCalEn;
0462         if (eeCalChi2 >= 0) {
0463           getCell(NNInputs::EcalChi2) += eeCalChi2 * eeCalEn;
0464           getCell(NNInputs::EcalEnergySumForPositiveChi2) += eeCalEn;
0465           getCell(NNInputs::EcalSizeForPositiveChi2) += 1.;
0466         }
0467       }
0468     }
0469 
0470     // calorechit_EB
0471     for (const auto& caloRecHit_eb : *caloRecHits.eb) {
0472       if (caloRecHit_eb.energy() <= 0)
0473         continue;
0474       const auto& position = caloRecHits.geometry->getGeometry(caloRecHit_eb.id())->getPosition();
0475       const float ebCalEn = caloRecHit_eb.energy();
0476       const float ebCalChi2 = caloRecHit_eb.chi2();
0477       if (reco::deltaR2(position, allTaus[tau_idx]->polarP4()) < dR2_max) {
0478         std::tie(deta, dphi, eta_idx, phi_idx) = getEtaPhiIndices(position, allTaus[tau_idx]->polarP4());
0479         getCell(NNInputs::EcalEnergySum) += ebCalEn;
0480         getCell(NNInputs::EcalSize) += 1.;
0481         getCell(NNInputs::EcalEnergyStdDev) += ebCalEn * ebCalEn;
0482         getCell(NNInputs::EcalDeltaEta) += deta * ebCalEn;
0483         getCell(NNInputs::EcalDeltaPhi) += dphi * ebCalEn;
0484         if (ebCalChi2 >= 0) {
0485           getCell(NNInputs::EcalChi2) += ebCalChi2 * ebCalEn;
0486           getCell(NNInputs::EcalEnergySumForPositiveChi2) += ebCalEn;
0487           getCell(NNInputs::EcalSizeForPositiveChi2) += 1.;
0488         }
0489       }
0490     }
0491 
0492     // calorechit_HBHE
0493     for (const auto& caloRecHit_hbhe : *caloRecHits.hbhe) {
0494       if (caloRecHit_hbhe.energy() <= 0)
0495         continue;
0496       const auto& position = caloRecHits.geometry->getGeometry(caloRecHit_hbhe.id())->getPosition();
0497       const float hbheCalEn = caloRecHit_hbhe.energy();
0498       const float hbheCalChi2 = caloRecHit_hbhe.chi2();
0499       if (reco::deltaR2(position, allTaus[tau_idx]->polarP4()) < dR2_max) {
0500         std::tie(deta, dphi, eta_idx, phi_idx) = getEtaPhiIndices(position, allTaus[tau_idx]->polarP4());
0501         getCell(NNInputs::HcalEnergySum) += hbheCalEn;
0502         getCell(NNInputs::HcalEnergyStdDev) += hbheCalEn * hbheCalEn;
0503         getCell(NNInputs::HcalSize) += 1.;
0504         getCell(NNInputs::HcalDeltaEta) += deta * hbheCalEn;
0505         getCell(NNInputs::HcalDeltaPhi) += dphi * hbheCalEn;
0506         if (hbheCalChi2 >= 0) {
0507           getCell(NNInputs::HcalChi2) += hbheCalChi2 * hbheCalEn;
0508           getCell(NNInputs::HcalEnergySumForPositiveChi2) += hbheCalEn;
0509           getCell(NNInputs::HcalSizeForPositiveChi2) += 1.;
0510         }
0511       }
0512     }
0513 
0514     // calorechit_HO
0515     for (const auto& caloRecHit_ho : *caloRecHits.ho) {
0516       if (caloRecHit_ho.energy() <= 0)
0517         continue;
0518       const auto& position = caloRecHits.geometry->getGeometry(caloRecHit_ho.id())->getPosition();
0519       const float hoCalEn = caloRecHit_ho.energy();
0520       if (reco::deltaR2(position, allTaus[tau_idx]->polarP4()) < dR2_max) {
0521         std::tie(deta, dphi, eta_idx, phi_idx) = getEtaPhiIndices(position, allTaus[tau_idx]->polarP4());
0522         getCell(NNInputs::HcalEnergySum) += hoCalEn;
0523         getCell(NNInputs::HcalEnergyStdDev) += hoCalEn * hoCalEn;
0524         getCell(NNInputs::HcalSize) += 1.;
0525         getCell(NNInputs::HcalDeltaEta) += deta * hoCalEn;
0526         getCell(NNInputs::HcalDeltaPhi) += dphi * hoCalEn;
0527       }
0528     }
0529 
0530     // normalize to sum and define stdDev
0531     for (eta_idx = 0; eta_idx < L2TauTagNNv1::nCellEta; eta_idx++) {
0532       for (phi_idx = 0; phi_idx < L2TauTagNNv1::nCellPhi; phi_idx++) {
0533         /* normalize eCal vars*/
0534         if (getCell(NNInputs::EcalEnergySum) > 0.) {
0535           getCell(NNInputs::EcalDeltaEta) /= getCell(NNInputs::EcalEnergySum);
0536           getCell(NNInputs::EcalDeltaPhi) /= getCell(NNInputs::EcalEnergySum);
0537         }
0538         if (getCell(NNInputs::EcalEnergySumForPositiveChi2) > 0.) {
0539           getCell(NNInputs::EcalChi2) /= getCell(NNInputs::EcalEnergySumForPositiveChi2);
0540         }
0541         if (getCell(NNInputs::EcalSize) > 1.) {
0542           // (stdDev - (enSum*enSum)/size) / (size-1)
0543           getCell(NNInputs::EcalEnergyStdDev) =
0544               (getCell(NNInputs::EcalEnergyStdDev) -
0545                (getCell(NNInputs::EcalEnergySum) * getCell(NNInputs::EcalEnergySum)) / getCell(NNInputs::EcalSize)) /
0546               (getCell(NNInputs::EcalSize) - 1);
0547         } else {
0548           getCell(NNInputs::EcalEnergyStdDev) = 0.;
0549         }
0550         /* normalize hCal Vars */
0551         if (getCell(NNInputs::HcalEnergySum) > 0.) {
0552           getCell(NNInputs::HcalDeltaEta) /= getCell(NNInputs::HcalEnergySum);
0553           getCell(NNInputs::HcalDeltaPhi) /= getCell(NNInputs::HcalEnergySum);
0554         }
0555         if (getCell(NNInputs::HcalEnergySumForPositiveChi2) > 0.) {
0556           getCell(NNInputs::HcalChi2) /= getCell(NNInputs::HcalEnergySumForPositiveChi2);
0557         }
0558         if (getCell(NNInputs::HcalSize) > 1.) {
0559           // (stdDev - (enSum*enSum)/size) / (size-1)
0560           getCell(NNInputs::HcalEnergyStdDev) =
0561               (getCell(NNInputs::HcalEnergyStdDev) -
0562                (getCell(NNInputs::HcalEnergySum) * getCell(NNInputs::HcalEnergySum)) / getCell(NNInputs::HcalSize)) /
0563               (getCell(NNInputs::HcalSize) - 1);
0564         } else {
0565           getCell(NNInputs::HcalEnergyStdDev) = 0.;
0566         }
0567       }
0568     }
0569   }
0570 }
0571 
0572 void L2TauNNProducer::selectGoodTracksAndVertices(const ZVertexSoA& patavtx_soa,
0573                                                   const pixelTrack::TrackSoA& patatracks_tsoa,
0574                                                   std::vector<int>& trkGood,
0575                                                   std::vector<int>& vtxGood) {
0576   const auto maxTracks = patatracks_tsoa.stride();
0577   const int nv = patavtx_soa.nvFinal;
0578   trkGood.clear();
0579   trkGood.reserve(maxTracks);
0580   vtxGood.clear();
0581   vtxGood.reserve(nv);
0582   auto const* quality = patatracks_tsoa.qualityData();
0583 
0584   // No need to sort either as the algorithms is just using the max (not even the location, just the max value of pt2sum).
0585   std::vector<float> pTSquaredSum(nv, 0);
0586   std::vector<int> nTrkAssociated(nv, 0);
0587 
0588   for (int32_t trk_idx = 0; trk_idx < maxTracks; ++trk_idx) {
0589     auto nHits = patatracks_tsoa.nHits(trk_idx);
0590     if (nHits == 0) {
0591       break;
0592     }
0593     int vtx_ass_to_track = patavtx_soa.idv[trk_idx];
0594     if (vtx_ass_to_track >= 0 && vtx_ass_to_track < nv) {
0595       auto patatrackPt = patatracks_tsoa.pt[trk_idx];
0596       ++nTrkAssociated[vtx_ass_to_track];
0597       if (patatrackPt >= trackPtMin_ && patatracks_tsoa.chi2(trk_idx) <= trackChi2Max_) {
0598         patatrackPt = std::min(patatrackPt, trackPtMax_);
0599         pTSquaredSum[vtx_ass_to_track] += patatrackPt * patatrackPt;
0600       }
0601     }
0602     if (nHits > 0 and quality[trk_idx] >= pixelTrack::Quality::loose) {
0603       trkGood.push_back(trk_idx);
0604     }
0605   }
0606   if (nv > 0) {
0607     const auto minFOM_fromFrac = (*std::max_element(pTSquaredSum.begin(), pTSquaredSum.end())) * fractionSumPt2_;
0608     for (int j = nv - 1; j >= 0 && vtxGood.size() < maxVtx_; --j) {
0609       auto vtx_idx = patavtx_soa.sortInd[j];
0610       assert(vtx_idx < nv);
0611       if (nTrkAssociated[vtx_idx] >= 2 && pTSquaredSum[vtx_idx] >= minFOM_fromFrac &&
0612           pTSquaredSum[vtx_idx] > minSumPt2_) {
0613         vtxGood.push_back(vtx_idx);
0614       }
0615     }
0616   }
0617 }
0618 
0619 std::pair<float, float> L2TauNNProducer::impactParameter(int it,
0620                                                          const pixelTrack::TrackSoA& patatracks_tsoa,
0621                                                          float patatrackPhi,
0622                                                          const reco::BeamSpot& beamspot,
0623                                                          const MagneticField* magfi) {
0624   auto const& fit = patatracks_tsoa.stateAtBS;
0625   /* dxy and dz */
0626   riemannFit::Vector5d ipar, opar;
0627   riemannFit::Matrix5d icov, ocov;
0628   fit.copyToDense(ipar, icov, it);
0629   riemannFit::transformToPerigeePlane(ipar, icov, opar, ocov);
0630   LocalTrajectoryParameters lpar(opar(0), opar(1), opar(2), opar(3), opar(4), 1.);
0631   float sp = std::sin(patatrackPhi);
0632   float cp = std::cos(patatrackPhi);
0633   Surface::RotationType Rotation(sp, -cp, 0, 0, 0, -1.f, cp, sp, 0);
0634   GlobalPoint BeamSpotPoint(beamspot.x0(), beamspot.y0(), beamspot.z0());
0635   Plane impPointPlane(BeamSpotPoint, Rotation);
0636   GlobalTrajectoryParameters gp(
0637       impPointPlane.toGlobal(lpar.position()), impPointPlane.toGlobal(lpar.momentum()), lpar.charge(), magfi);
0638   GlobalPoint vv = gp.position();
0639   math::XYZPoint pos(vv.x(), vv.y(), vv.z());
0640   GlobalVector pp = gp.momentum();
0641   math::XYZVector mom(pp.x(), pp.y(), pp.z());
0642   auto lambda = M_PI_2 - pp.theta();
0643   auto phi = pp.phi();
0644   float patatrackDxy = -vv.x() * std::sin(phi) + vv.y() * std::cos(phi);
0645   float patatrackDz =
0646       (vv.z() * std::cos(lambda) - (vv.x() * std::cos(phi) + vv.y() * std::sin(phi)) * std::sin(lambda)) /
0647       std::cos(lambda);
0648   return std::make_pair(patatrackDxy, patatrackDz);
0649 }
0650 
0651 void L2TauNNProducer::fillPatatracks(tensorflow::Tensor& cellGridMatrix,
0652                                      const std::vector<l1t::TauRef>& allTaus,
0653                                      const pixelTrack::TrackSoA& patatracks_tsoa,
0654                                      const ZVertexSoA& patavtx_soa,
0655                                      const reco::BeamSpot& beamspot,
0656                                      const MagneticField* magfi) {
0657   using NNInputs = L2TauTagNNv1::NNInputs;
0658   float deta, dphi;
0659   int eta_idx = 0;
0660   int phi_idx = 0;
0661   int tau_idx = 0;
0662 
0663   auto getCell = [&](NNInputs input) -> float& {
0664     return getCellImpl(cellGridMatrix, tau_idx, phi_idx, eta_idx, input);
0665   };
0666 
0667   std::vector<int> trkGood;
0668   std::vector<int> vtxGood;
0669 
0670   selectGoodTracksAndVertices(patavtx_soa, patatracks_tsoa, trkGood, vtxGood);
0671 
0672   const int nTaus = allTaus.size();
0673   for (tau_idx = 0; tau_idx < nTaus; tau_idx++) {
0674     const float tauEta = allTaus[tau_idx]->eta();
0675     const float tauPhi = allTaus[tau_idx]->phi();
0676 
0677     for (const auto it : trkGood) {
0678       const float patatrackPt = patatracks_tsoa.pt[it];
0679       if (patatrackPt <= 0)
0680         continue;
0681       const float patatrackPhi = patatracks_tsoa.phi(it);
0682       const float patatrackEta = patatracks_tsoa.eta(it);
0683       const float patatrackCharge = patatracks_tsoa.charge(it);
0684       const float patatrackChi2OverNdof = patatracks_tsoa.chi2(it);
0685       const auto nHits = patatracks_tsoa.nHits(it);
0686       if (nHits <= 0)
0687         continue;
0688       const int patatrackNdof = 2 * std::min(6, nHits) - 5;
0689 
0690       const int vtx_idx_assTrk = patavtx_soa.idv[it];
0691       if (reco::deltaR2(patatrackEta, patatrackPhi, tauEta, tauPhi) < dR2_max) {
0692         std::tie(deta, dphi, eta_idx, phi_idx) =
0693             getEtaPhiIndices(patatrackEta, patatrackPhi, allTaus[tau_idx]->polarP4());
0694         getCell(NNInputs::PatatrackPtSum) += patatrackPt;
0695         getCell(NNInputs::PatatrackSize) += 1.;
0696         getCell(NNInputs::PatatrackChargeSum) += patatrackCharge;
0697         getCell(NNInputs::PatatrackDeltaEta) += deta * patatrackPt;
0698         getCell(NNInputs::PatatrackDeltaPhi) += dphi * patatrackPt;
0699         getCell(NNInputs::PatatrackChi2OverNdof) += patatrackChi2OverNdof * patatrackPt;
0700         getCell(NNInputs::PatatrackNdof) += patatrackNdof * patatrackPt;
0701         std::pair<float, float> impactParameters = impactParameter(it, patatracks_tsoa, patatrackPhi, beamspot, magfi);
0702         getCell(NNInputs::PatatrackDxy) += impactParameters.first * patatrackPt;
0703         getCell(NNInputs::PatatrackDz) += impactParameters.second * patatrackPt;
0704         if ((std::find(vtxGood.begin(), vtxGood.end(), vtx_idx_assTrk) != vtxGood.end())) {
0705           getCell(NNInputs::PatatrackPtSumWithVertex) += patatrackPt;
0706           getCell(NNInputs::PatatrackSizeWithVertex) += 1.;
0707         }
0708       }
0709     }
0710 
0711     // normalize to sum and define stdDev
0712     for (eta_idx = 0; eta_idx < L2TauTagNNv1::nCellEta; eta_idx++) {
0713       for (phi_idx = 0; phi_idx < L2TauTagNNv1::nCellPhi; phi_idx++) {
0714         getCell(NNInputs::nVertices) = vtxGood.size();
0715         if (getCell(NNInputs::PatatrackPtSum) > 0.) {
0716           getCell(NNInputs::PatatrackDeltaEta) /= getCell(NNInputs::PatatrackPtSum);
0717           getCell(NNInputs::PatatrackDeltaPhi) /= getCell(NNInputs::PatatrackPtSum);
0718           getCell(NNInputs::PatatrackChi2OverNdof) /= getCell(NNInputs::PatatrackPtSum);
0719           getCell(NNInputs::PatatrackNdof) /= getCell(NNInputs::PatatrackPtSum);
0720           getCell(NNInputs::PatatrackDxy) /= getCell(NNInputs::PatatrackPtSum);
0721           getCell(NNInputs::PatatrackDz) /= getCell(NNInputs::PatatrackPtSum);
0722         }
0723       }
0724     }
0725   }
0726 }
0727 
0728 std::vector<float> L2TauNNProducer::getTauScore(const tensorflow::Tensor& cellGridMatrix) {
0729   std::vector<tensorflow::Tensor> pred_tensor;
0730   tensorflow::run(L2cacheData_->session, {{inputTensorName_, cellGridMatrix}}, {outputTensorName_}, &pred_tensor);
0731   const int nTau = cellGridMatrix.shape().dim_size(0);
0732   std::vector<float> pred_vector(nTau);
0733   for (int tau_idx = 0; tau_idx < nTau; ++tau_idx) {
0734     pred_vector[tau_idx] = pred_tensor[0].matrix<float>()(tau_idx, 0);
0735   }
0736 
0737   return pred_vector;
0738 }
0739 
0740 void L2TauNNProducer::produce(edm::Event& event, const edm::EventSetup& eventsetup) {
0741   std::vector<std::vector<size_t>> TauCollectionMap(L1TauDesc_.size());
0742   l1t::TauVectorRef allTaus;
0743 
0744   for (size_t inp_idx = 0; inp_idx < L1TauDesc_.size(); inp_idx++) {
0745     l1t::TauVectorRef l1Taus;
0746     auto const& l1TriggeredTaus = event.get(L1TauDesc_[inp_idx].inputToken_);
0747     l1TriggeredTaus.getObjects(trigger::TriggerL1Tau, l1Taus);
0748     TauCollectionMap.at(inp_idx).resize(l1Taus.size());
0749 
0750     for (size_t l1_idx = 0; l1_idx < l1Taus.size(); l1_idx++) {
0751       size_t tau_idx;
0752       const auto iter = std::find(allTaus.begin(), allTaus.end(), l1Taus[l1_idx]);
0753       if (iter != allTaus.end()) {
0754         tau_idx = std::distance(allTaus.begin(), iter);
0755       } else {
0756         allTaus.push_back(l1Taus[l1_idx]);
0757         tau_idx = allTaus.size() - 1;
0758       }
0759       TauCollectionMap.at(inp_idx).at(l1_idx) = tau_idx;
0760     }
0761   }
0762   const auto ebCal = event.getHandle(ebToken_);
0763   const auto eeCal = event.getHandle(eeToken_);
0764   const auto hbhe = event.getHandle(hbheToken_);
0765   const auto ho = event.getHandle(hoToken_);
0766   const auto& patatracks_SoA = *event.get(pataTracksToken_);
0767   const auto& vertices_SoA = *event.get(pataVerticesToken_);
0768   const auto bsHandle = event.getHandle(beamSpotToken_);
0769 
0770   auto const fieldESH = eventsetup.getHandle(bFieldToken_);
0771   auto const geometry = eventsetup.getHandle(geometryToken_);
0772 
0773   caloRecHitCollections caloRecHits;
0774   caloRecHits.hbhe = &*hbhe;
0775   caloRecHits.ho = &*ho;
0776   caloRecHits.eb = &*ebCal;
0777   caloRecHits.ee = &*eeCal;
0778   caloRecHits.geometry = &*geometry;
0779 
0780   const int nTaus = allTaus.size();
0781   tensorflow::Tensor cellGridMatrix(tensorflow::DT_FLOAT,
0782                                     {nTaus, L2TauTagNNv1::nCellEta, L2TauTagNNv1::nCellPhi, L2TauTagNNv1::nVars});
0783   const int n_inputs = nTaus * L2TauTagNNv1::nCellEta * L2TauTagNNv1::nCellPhi * L2TauTagNNv1::nVars;
0784   for (int input_idx = 0; input_idx < n_inputs; ++input_idx) {
0785     cellGridMatrix.flat<float>()(input_idx) = 0;
0786   }
0787   fillL1TauVars(cellGridMatrix, allTaus);
0788 
0789   fillCaloRecHits(cellGridMatrix, allTaus, caloRecHits);
0790 
0791   fillPatatracks(cellGridMatrix, allTaus, patatracks_SoA, vertices_SoA, *bsHandle, fieldESH.product());
0792 
0793   standardizeTensor(cellGridMatrix);
0794 
0795   if (debugLevel_ > 0) {
0796     checknan(cellGridMatrix, debugLevel_);
0797   }
0798 
0799   std::vector<float> tau_score = getTauScore(cellGridMatrix);
0800 
0801   for (size_t inp_idx = 0; inp_idx < L1TauDesc_.size(); inp_idx++) {
0802     const size_t nTau = TauCollectionMap[inp_idx].size();
0803     auto tau_tags = std::make_unique<std::vector<float>>(nTau);
0804     for (size_t tau_pos = 0; tau_pos < nTau; ++tau_pos) {
0805       const auto tau_idx = TauCollectionMap[inp_idx][tau_pos];
0806       if (debugLevel_ > 0) {
0807         edm::LogInfo("DebugInfo") << event.id().event() << " \t " << (allTaus[tau_idx])->pt() << " \t "
0808                                   << tau_score.at(tau_idx) << std::endl;
0809       }
0810       (*tau_tags)[tau_pos] = tau_score.at(tau_idx);
0811     }
0812     event.put(std::move(tau_tags), L1TauDesc_[inp_idx].CollectionName);
0813   }
0814 }
0815 //define this as a plug-in
0816 #include "FWCore/Framework/interface/MakerMacros.h"
0817 DEFINE_FWK_MODULE(L2TauNNProducer);