Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:36:49 

0001 // -*- C++ -*-
0002 //
0003 // Package: L1CaloTrigger
0004 // Class: L1CaloJetProducer
0005 //
0006 /**\class L1CaloJetProducer L1CaloJetProducer.cc
0007 
0008 Description: 
0009 Beginning with HCAL TPs, create HCAL jet, then
0010 take L1EG crystal clusters from L1EGammaCrystalsProducer.cc
0011 and clusters them within fixed number of trigger towers
0012 
0013 Implementation:
0014 [Notes on implementation]
0015 */
0016 //
0017 // Original Author: Tyler Ruggles
0018 // Created: Tue Aug 29 2018
0019 // $Id$
0020 //
0021 //
0022 
0023 #include "DataFormats/Math/interface/deltaR.h"
0024 #include "DataFormats/Math/interface/deltaPhi.h"
0025 
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/stream/EDProducer.h"
0028 #include "FWCore/Framework/interface/ESHandle.h"
0029 #include "FWCore/ServiceRegistry/interface/Service.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0033 
0034 #include <iostream>
0035 
0036 #include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h"
0037 //#include "DataFormats/L1TCalorimeterPhase2/interface/CaloCrystalCluster.h"
0038 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloJet.h"
0039 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloTower.h"
0040 #include "DataFormats/L1THGCal/interface/HGCalTower.h"
0041 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0042 
0043 #include "CalibFormats/CaloTPG/interface/CaloTPGTranscoder.h"
0044 #include "CalibFormats/CaloTPG/interface/CaloTPGRecord.h"
0045 
0046 #include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
0047 
0048 // For pT calibrations
0049 #include "TF1.h"
0050 
0051 // Run2/PhaseI output formats
0052 #include "DataFormats/L1Trigger/interface/Tau.h"
0053 #include "DataFormats/L1Trigger/interface/Jet.h"
0054 
0055 class L1CaloJetProducer : public edm::stream::EDProducer<> {
0056 public:
0057   explicit L1CaloJetProducer(const edm::ParameterSet &);
0058 
0059 private:
0060   void produce(edm::Event &, const edm::EventSetup &) override;
0061   //bool cluster_passes_base_cuts(float &cluster_pt, float &cluster_eta, float &iso, float &e2x5, float &e5x5) const;
0062   int tower_diPhi(int &iPhi_1, int &iPhi_2) const;
0063   int tower_diEta(int &iEta_1, int &iEta_2) const;
0064   float get_deltaR(reco::Candidate::PolarLorentzVector &p4_1, reco::Candidate::PolarLorentzVector &p4_2) const;
0065   float get_hcal_calibration(float &jet_pt, float &ecal_pt, float &ecal_L1EG_jet_pt, float &jet_eta) const;
0066   float get_tau_pt_calibration(float &tau_pt, float &ecal_pt, float &l1EG_pt, float &n_L1EGs, float &tau_eta) const;
0067   int loose_iso_tau_wp(float &tau_pt, float &tau_iso_et, float &tau_eta) const;
0068 
0069   double HcalTpEtMin;
0070   double EcalTpEtMin;
0071   double HGCalHadTpEtMin;
0072   double HGCalEmTpEtMin;
0073   double HFTpEtMin;
0074   double EtMinForSeedHit;
0075   double EtMinForCollection;
0076   double EtMinForTauCollection;
0077 
0078   // For fetching jet calibrations
0079   std::vector<double> jetPtBins;
0080   std::vector<double> emFractionBinsBarrel;
0081   std::vector<double> absEtaBinsBarrel;
0082   std::vector<double> jetCalibrationsBarrel;
0083   std::vector<double> emFractionBinsHGCal;
0084   std::vector<double> absEtaBinsHGCal;
0085   std::vector<double> jetCalibrationsHGCal;
0086   std::vector<double> emFractionBinsHF;
0087   std::vector<double> absEtaBinsHF;
0088   std::vector<double> jetCalibrationsHF;
0089 
0090   // For fetching tau calibrations
0091   std::vector<double> tauPtBins;
0092   std::vector<double> tauAbsEtaBinsBarrel;
0093   std::vector<double> tauCalibrationsBarrel;
0094   std::vector<edm::ParameterSet> tauL1egInfoBarrel;
0095   std::vector<double> tauAbsEtaBinsHGCal;
0096   std::vector<double> tauCalibrationsHGCal;
0097   std::vector<edm::ParameterSet> tauL1egInfoHGCal;
0098 
0099   // For storing jet calibrations
0100   std::vector<std::vector<std::vector<double>>> calibrationsBarrel;
0101   std::vector<std::vector<std::vector<double>>> calibrationsHGCal;
0102   std::vector<std::vector<std::vector<double>>> calibrationsHF;
0103 
0104   // For storing tau calibration info
0105   std::map<double, std::vector<double>> tauL1egInfoMapBarrel;
0106   std::map<double, std::vector<double>> tauL1egInfoMapHGCal;
0107   std::vector<double> tauL1egValuesBarrel;  // To preserve ordering
0108   std::vector<double> tauL1egValuesHGCal;   // To preserve ordering
0109   std::vector<std::vector<std::vector<std::vector<double>>>> tauPtCalibrationsBarrel;
0110   std::vector<std::vector<std::vector<std::vector<double>>>> tauPtCalibrationsHGCal;
0111 
0112   bool debug;
0113   edm::EDGetTokenT<l1tp2::CaloTowerCollection> l1TowerToken_;
0114   edm::Handle<l1tp2::CaloTowerCollection> l1CaloTowerHandle;
0115 
0116   // TF1s defining tau isolation thresholds
0117   TF1 isoTauBarrel = TF1("isoTauBarrelFunction", "([0] + [1]*TMath::Exp(-[2]*x))");
0118   TF1 isoTauHGCal = TF1("isoTauHGCalFunction", "([0] + [1]*TMath::Exp(-[2]*x))");
0119 
0120   class l1CaloJetObj {
0121   public:
0122     bool barrelSeeded = true;  // default to barrel seeded
0123     reco::Candidate::PolarLorentzVector jetCluster;
0124     reco::Candidate::PolarLorentzVector hcalJetCluster;
0125     reco::Candidate::PolarLorentzVector ecalJetCluster;
0126     reco::Candidate::PolarLorentzVector seedTower;
0127     //reco::Candidate::PolarLorentzVector leadingL1EG;
0128     reco::Candidate::PolarLorentzVector l1egJetCluster;
0129     float jetClusterET = 0.;
0130     float hcalJetClusterET = 0.;
0131     float ecalJetClusterET = 0.;
0132     float seedTowerET = 0.;
0133     float leadingL1EGET = 0.;
0134     float l1egJetClusterET = 0.;
0135 
0136     // For decay mode related checks with CaloTaus
0137     std::vector<std::vector<float>> associated_l1EGs_;
0138 
0139     int seed_iEta = -99;
0140     int seed_iPhi = -99;
0141 
0142     float hcal_seed = 0.;
0143     //float hcal_3x3 = 0.;
0144     float hcal_3x5 = 0.;
0145     //float hcal_5x5 = 0.;
0146     //float hcal_5x7 = 0.;
0147     float hcal_7x7 = 0.;
0148     //float hcal_2x3 = 0.;
0149     //float hcal_2x3_1 = 0.;
0150     //float hcal_2x3_2 = 0.;
0151     //float hcal_2x2 = 0.;
0152     //float hcal_2x2_1 = 0.;
0153     //float hcal_2x2_2 = 0.;
0154     //float hcal_2x2_3 = 0.;
0155     //float hcal_2x2_4 = 0.;
0156     float hcal_nHits = 0.;
0157 
0158     float ecal_seed = 0.;
0159     //float ecal_3x3 = 0.;
0160     float ecal_3x5 = 0.;
0161     //float ecal_5x5 = 0.;
0162     //float ecal_5x7 = 0.;
0163     float ecal_7x7 = 0.;
0164     //float ecal_2x3 = 0.;
0165     //float ecal_2x3_1 = 0.;
0166     //float ecal_2x3_2 = 0.;
0167     //float ecal_2x2 = 0.;
0168     //float ecal_2x2_1 = 0.;
0169     //float ecal_2x2_2 = 0.;
0170     //float ecal_2x2_3 = 0.;
0171     //float ecal_2x2_4 = 0.;
0172     float ecal_nHits = 0.;
0173 
0174     float l1eg_seed = 0.;
0175     //float l1eg_3x3 = 0.;
0176     float l1eg_3x5 = 0.;
0177     //float l1eg_5x5 = 0.;
0178     //float l1eg_5x7 = 0.;
0179     float l1eg_7x7 = 0.;
0180     //float l1eg_2x3 = 0.;
0181     //float l1eg_2x3_1 = 0.;
0182     //float l1eg_2x3_2 = 0.;
0183     //float l1eg_2x2 = 0.;
0184     //float l1eg_2x2_1 = 0.;
0185     //float l1eg_2x2_2 = 0.;
0186     //float l1eg_2x2_3 = 0.;
0187     //float l1eg_2x2_4 = 0.;
0188     float l1eg_nHits = 0.;
0189     float n_l1eg_HoverE_LessThreshold = 0.;
0190 
0191     float l1eg_nL1EGs = 0.;
0192     float l1eg_nL1EGs_standaloneSS = 0.;
0193     float l1eg_nL1EGs_standaloneIso = 0.;
0194     float l1eg_nL1EGs_trkMatchSS = 0.;
0195     float l1eg_nL1EGs_trkMatchIso = 0.;
0196 
0197     float total_seed = 0.;
0198     //float total_3x3 = 0.;
0199     float total_3x5 = 0.;
0200     //float total_5x5 = 0.;
0201     //float total_5x7 = 0.;
0202     float total_7x7 = 0.;
0203     //float total_2x3 = 0.;
0204     //float total_2x3_1 = 0.;
0205     //float total_2x3_2 = 0.;
0206     //float total_2x2 = 0.;
0207     //float total_2x2_1 = 0.;
0208     //float total_2x2_2 = 0.;
0209     //float total_2x2_3 = 0.;
0210     //float total_2x2_4 = 0.;
0211     float total_nHits = 0.;
0212 
0213     void Init() {
0214       SetJetClusterP4(0., 0., 0., 0.);
0215       SetHcalJetClusterP4(0., 0., 0., 0.);
0216       SetEcalJetClusterP4(0., 0., 0., 0.);
0217       SetSeedP4(0., 0., 0., 0.);
0218       //SetLeadingL1EGP4( 0., 0., 0., 0. );
0219       SetL1EGJetP4(0., 0., 0., 0.);
0220     }
0221 
0222     void SetJetClusterP4(double pt, double eta, double phi, double mass) {
0223       this->jetCluster.SetPt(pt);
0224       this->jetCluster.SetEta(eta);
0225       this->jetCluster.SetPhi(phi);
0226       this->jetCluster.SetM(mass);
0227     }
0228     void SetHcalJetClusterP4(double pt, double eta, double phi, double mass) {
0229       this->hcalJetCluster.SetPt(pt);
0230       this->hcalJetCluster.SetEta(eta);
0231       this->hcalJetCluster.SetPhi(phi);
0232       this->hcalJetCluster.SetM(mass);
0233     }
0234     void SetEcalJetClusterP4(double pt, double eta, double phi, double mass) {
0235       this->ecalJetCluster.SetPt(pt);
0236       this->ecalJetCluster.SetEta(eta);
0237       this->ecalJetCluster.SetPhi(phi);
0238       this->ecalJetCluster.SetM(mass);
0239     }
0240     void SetSeedP4(double pt, double eta, double phi, double mass) {
0241       this->seedTower.SetPt(pt);
0242       this->seedTower.SetEta(eta);
0243       this->seedTower.SetPhi(phi);
0244       this->seedTower.SetM(mass);
0245     }
0246     //void SetLeadingL1EGP4( double pt, double eta, double phi, double mass )
0247     //{
0248     //    this->leadingL1EG.SetPt( pt );
0249     //    this->leadingL1EG.SetEta( eta );
0250     //    this->leadingL1EG.SetPhi( phi );
0251     //    this->leadingL1EG.SetM( mass );
0252     //}
0253     void SetL1EGJetP4(double pt, double eta, double phi, double mass) {
0254       this->l1egJetCluster.SetPt(pt);
0255       this->l1egJetCluster.SetEta(eta);
0256       this->l1egJetCluster.SetPhi(phi);
0257       this->l1egJetCluster.SetM(mass);
0258     }
0259   };
0260 
0261   class simpleL1obj {
0262   public:
0263     bool stale = false;  // Hits become stale once used in clustering algorithm to prevent overlap in clusters
0264     bool associated_with_tower =
0265         false;                         // L1EGs become associated with a tower to find highest ET total for seeding jets
0266     bool passesStandaloneSS = false;   // Store whether any of the portions of a WP are passed
0267     bool passesStandaloneIso = false;  // Store whether any of the portions of a WP are passed
0268     bool passesTrkMatchSS = false;     // Store whether any of the portions of a WP are passed
0269     bool passesTrkMatchIso = false;    // Store whether any of the portions of a WP are passed
0270     reco::Candidate::PolarLorentzVector p4;
0271 
0272     void SetP4(double pt, double eta, double phi, double mass) {
0273       this->p4.SetPt(pt);
0274       this->p4.SetEta(eta);
0275       this->p4.SetPhi(phi);
0276       this->p4.SetM(mass);
0277     }
0278     inline float pt() const { return p4.pt(); };
0279     inline float eta() const { return p4.eta(); };
0280     inline float phi() const { return p4.phi(); };
0281     inline float M() const { return p4.M(); };
0282     inline reco::Candidate::PolarLorentzVector GetP4() const { return p4; };
0283   };
0284 
0285   class SimpleCaloHit {
0286   public:
0287     int towerIEta = -99;
0288     int towerIPhi = -99;
0289     float towerEta = -99;
0290     float towerPhi = -99;
0291     float ecalTowerEt = 0.;
0292     float hcalTowerEt = 0.;
0293     float l1egTowerEt = 0.;
0294     float total_tower_et = 0.;
0295     //float total_tower_plus_L1EGs_et=0.;
0296     bool stale = false;    // Hits become stale once used in clustering algorithm to prevent overlap in clusters
0297     bool isBarrel = true;  // Defaults to a barrel hit
0298 
0299     // L1EG info
0300     int nL1eg = 0;
0301     int l1egTrkSS = 0;
0302     int l1egTrkIso = 0;
0303     int l1egStandaloneSS = 0;
0304     int l1egStandaloneIso = 0;
0305   };
0306 };
0307 
0308 L1CaloJetProducer::L1CaloJetProducer(const edm::ParameterSet &iConfig)
0309     : HcalTpEtMin(iConfig.getParameter<double>("HcalTpEtMin")),                      // Should default to 0 MeV
0310       EcalTpEtMin(iConfig.getParameter<double>("EcalTpEtMin")),                      // Should default to 0 MeV
0311       HGCalHadTpEtMin(iConfig.getParameter<double>("HGCalHadTpEtMin")),              // Should default to 0 MeV
0312       HGCalEmTpEtMin(iConfig.getParameter<double>("HGCalEmTpEtMin")),                // Should default to 0 MeV
0313       HFTpEtMin(iConfig.getParameter<double>("HFTpEtMin")),                          // Should default to 0 MeV
0314       EtMinForSeedHit(iConfig.getParameter<double>("EtMinForSeedHit")),              // Should default to 2.5 GeV
0315       EtMinForCollection(iConfig.getParameter<double>("EtMinForCollection")),        // Testing 10 GeV
0316       EtMinForTauCollection(iConfig.getParameter<double>("EtMinForTauCollection")),  // Testing 10 GeV
0317       jetPtBins(iConfig.getParameter<std::vector<double>>("jetPtBins")),
0318       emFractionBinsBarrel(iConfig.getParameter<std::vector<double>>("emFractionBinsBarrel")),
0319       absEtaBinsBarrel(iConfig.getParameter<std::vector<double>>("absEtaBinsBarrel")),
0320       jetCalibrationsBarrel(iConfig.getParameter<std::vector<double>>("jetCalibrationsBarrel")),
0321       emFractionBinsHGCal(iConfig.getParameter<std::vector<double>>("emFractionBinsHGCal")),
0322       absEtaBinsHGCal(iConfig.getParameter<std::vector<double>>("absEtaBinsHGCal")),
0323       jetCalibrationsHGCal(iConfig.getParameter<std::vector<double>>("jetCalibrationsHGCal")),
0324       emFractionBinsHF(iConfig.getParameter<std::vector<double>>("emFractionBinsHF")),
0325       absEtaBinsHF(iConfig.getParameter<std::vector<double>>("absEtaBinsHF")),
0326       jetCalibrationsHF(iConfig.getParameter<std::vector<double>>("jetCalibrationsHF")),
0327       tauPtBins(iConfig.getParameter<std::vector<double>>("tauPtBins")),
0328       tauAbsEtaBinsBarrel(iConfig.getParameter<std::vector<double>>("tauAbsEtaBinsBarrel")),
0329       tauCalibrationsBarrel(iConfig.getParameter<std::vector<double>>("tauCalibrationsBarrel")),
0330       tauL1egInfoBarrel(iConfig.getParameter<std::vector<edm::ParameterSet>>("tauL1egInfoBarrel")),
0331       tauAbsEtaBinsHGCal(iConfig.getParameter<std::vector<double>>("tauAbsEtaBinsHGCal")),
0332       tauCalibrationsHGCal(iConfig.getParameter<std::vector<double>>("tauCalibrationsHGCal")),
0333       tauL1egInfoHGCal(iConfig.getParameter<std::vector<edm::ParameterSet>>("tauL1egInfoHGCal")),
0334       debug(iConfig.getParameter<bool>("debug")),
0335       l1TowerToken_(consumes<l1tp2::CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("l1CaloTowers")))
0336 
0337 {
0338   produces<l1tp2::CaloJetsCollection>("L1CaloJetsNoCuts");
0339   //produces<l1tp2::CaloJetsCollection>("L1CaloJetsWithCuts");
0340   //produces<l1extra::L1JetParticleCollection>("L1CaloClusterCollectionWithCuts");
0341   produces<BXVector<l1t::Jet>>("L1CaloJetCollectionBXV");
0342   produces<BXVector<l1t::Tau>>("L1CaloTauCollectionBXV");
0343 
0344   if (debug) {
0345     LogDebug("L1CaloJetProducer") << " HcalTpEtMin = " << HcalTpEtMin << " EcalTpEtMin = " << EcalTpEtMin << "\n";
0346   }
0347 
0348   // Fill the calibration 3D vector
0349   // Dimension 1 is AbsEta bin
0350   // Dimension 2 is EM Fraction bin
0351   // Dimension 3 is jet pT bin which is filled with the actual callibration value
0352   // size()-1 b/c the inputs have lower and upper bounds
0353   // Do Barrel, then HGCal, then HF
0354   int index = 0;
0355   //calibrations[em_frac][abs_eta].push_back( jetCalibrationsBarrel.at(index) );
0356   for (unsigned int abs_eta = 0; abs_eta < absEtaBinsBarrel.size() - 1; abs_eta++) {
0357     std::vector<std::vector<double>> em_bins;
0358     for (unsigned int em_frac = 0; em_frac < emFractionBinsBarrel.size() - 1; em_frac++) {
0359       std::vector<double> pt_bin_calibs;
0360       for (unsigned int pt = 0; pt < jetPtBins.size() - 1; pt++) {
0361         pt_bin_calibs.push_back(jetCalibrationsBarrel.at(index));
0362         index++;
0363       }
0364       em_bins.push_back(pt_bin_calibs);
0365     }
0366     calibrationsBarrel.push_back(em_bins);
0367   }
0368   if (debug) {
0369     LogDebug("L1CaloJetProducer") << " Loading Barrel calibrations: Loaded " << index
0370                                   << " values vs. size() of input calibration file: "
0371                                   << int(jetCalibrationsBarrel.size()) << "\n";
0372   }
0373 
0374   index = 0;
0375   //calibrations[em_frac][abs_eta].push_back( jetCalibrationsHGCal.at(index) );
0376   for (unsigned int abs_eta = 0; abs_eta < absEtaBinsHGCal.size() - 1; abs_eta++) {
0377     std::vector<std::vector<double>> em_bins;
0378     for (unsigned int em_frac = 0; em_frac < emFractionBinsHGCal.size() - 1; em_frac++) {
0379       std::vector<double> pt_bin_calibs;
0380       for (unsigned int pt = 0; pt < jetPtBins.size() - 1; pt++) {
0381         pt_bin_calibs.push_back(jetCalibrationsHGCal.at(index));
0382         index++;
0383       }
0384       em_bins.push_back(pt_bin_calibs);
0385     }
0386     calibrationsHGCal.push_back(em_bins);
0387   }
0388   if (debug) {
0389     LogDebug("L1CaloJetProducer") << " Loading HGCal calibrations: Loaded " << index
0390                                   << " values vs. size() of input calibration file: "
0391                                   << int(jetCalibrationsHGCal.size()) << "\n";
0392   }
0393 
0394   index = 0;
0395   //calibrations[em_frac][abs_eta].push_back( jetCalibrationsHF.at(index) );
0396   for (unsigned int abs_eta = 0; abs_eta < absEtaBinsHF.size() - 1; abs_eta++) {
0397     std::vector<std::vector<double>> em_bins;
0398     for (unsigned int em_frac = 0; em_frac < emFractionBinsHF.size() - 1; em_frac++) {
0399       std::vector<double> pt_bin_calibs;
0400       for (unsigned int pt = 0; pt < jetPtBins.size() - 1; pt++) {
0401         pt_bin_calibs.push_back(jetCalibrationsHF.at(index));
0402         index++;
0403       }
0404       em_bins.push_back(pt_bin_calibs);
0405     }
0406     calibrationsHF.push_back(em_bins);
0407   }
0408   if (debug) {
0409     LogDebug("L1CaloJetProducer") << " Loading HF calibrations: Loaded " << index
0410                                   << " values vs. size() of input calibration file: " << int(jetCalibrationsHF.size())
0411                                   << "\n";
0412   }
0413 
0414   // Load Tau L1EG-base calibration info into maps
0415   for (auto &first : tauL1egInfoBarrel) {
0416     if (debug) {
0417       LogDebug("L1CaloJetProducer") << " barrel l1egCount = " << first.getParameter<double>("l1egCount") << "\n";
0418       for (auto &em_frac : first.getParameter<std::vector<double>>("l1egEmFractions")) {
0419         LogDebug("L1CaloJetProducer") << " - EM = " << em_frac << "\n";
0420       }
0421     }
0422     double l1egCount = first.getParameter<double>("l1egCount");
0423     std::vector<double> l1egEmFractions = first.getParameter<std::vector<double>>("l1egEmFractions");
0424     tauL1egInfoMapBarrel[l1egCount] = l1egEmFractions;
0425     tauL1egValuesBarrel.push_back(l1egCount);
0426   }
0427   std::sort(tauL1egValuesBarrel.begin(), tauL1egValuesBarrel.end());
0428   for (auto &first : tauL1egInfoHGCal) {
0429     if (debug) {
0430       LogDebug("L1CaloJetProducer") << " hgcal l1egCount = " << first.getParameter<double>("l1egCount") << "\n";
0431       for (auto &em_frac : first.getParameter<std::vector<double>>("l1egEmFractions")) {
0432         LogDebug("L1CaloJetProducer") << " - EM = " << em_frac << "\n";
0433       }
0434     }
0435     double l1egCount = first.getParameter<double>("l1egCount");
0436     std::vector<double> l1egEmFractions = first.getParameter<std::vector<double>>("l1egEmFractions");
0437     tauL1egInfoMapHGCal[l1egCount] = l1egEmFractions;
0438     tauL1egValuesHGCal.push_back(l1egCount);
0439   }
0440   std::sort(tauL1egValuesHGCal.begin(), tauL1egValuesHGCal.end());
0441   // Fill the calibration 4D vector
0442   // Dimension 1 is AbsEta bin
0443   // Dimension 2 is L1EG count
0444   // Dimension 3 is EM Fraction bin
0445   // Dimension 4 is tau pT bin which is filled with the actual callibration value
0446   // size()-1 b/c the inputs have lower and upper bounds (except L1EG b/c that is a cound)
0447   // Do Barrel, then HGCal
0448   //
0449   // Note to future developers: be very concious of the order in which the calibrations are printed
0450   // out in tool which makse the cfg files.  You need to match that exactly when loading them and
0451   // using the calibrations below.
0452   index = 0;
0453   for (unsigned int abs_eta = 0; abs_eta < tauAbsEtaBinsBarrel.size() - 1; abs_eta++) {
0454     std::vector<std::vector<std::vector<double>>> l1eg_bins;
0455     for (auto &l1eg_info : tauL1egInfoMapBarrel) {
0456       std::vector<std::vector<double>> em_bins;
0457       for (unsigned int em_frac = 0; em_frac < l1eg_info.second.size() - 1; em_frac++) {
0458         std::vector<double> pt_bin_calibs;
0459         for (unsigned int pt = 0; pt < tauPtBins.size() - 1; pt++) {
0460           pt_bin_calibs.push_back(tauCalibrationsBarrel.at(index));
0461           index++;
0462         }
0463         em_bins.push_back(pt_bin_calibs);
0464       }
0465       l1eg_bins.push_back(em_bins);
0466     }
0467     tauPtCalibrationsBarrel.push_back(l1eg_bins);
0468   }
0469   if (debug) {
0470     LogDebug("L1CaloJetProducer") << " Loading Barrel calibrations: Loaded " << index
0471                                   << " values vs. size() of input calibration file: "
0472                                   << int(tauCalibrationsBarrel.size()) << "\n";
0473   }
0474 
0475   index = 0;
0476   for (unsigned int abs_eta = 0; abs_eta < tauAbsEtaBinsHGCal.size() - 1; abs_eta++) {
0477     std::vector<std::vector<std::vector<double>>> l1eg_bins;
0478     for (const auto &l1eg_info : tauL1egInfoMapHGCal) {
0479       std::vector<std::vector<double>> em_bins;
0480       for (unsigned int em_frac = 0; em_frac < l1eg_info.second.size() - 1; em_frac++) {
0481         std::vector<double> pt_bin_calibs;
0482         for (unsigned int pt = 0; pt < tauPtBins.size() - 1; pt++) {
0483           pt_bin_calibs.push_back(tauCalibrationsHGCal.at(index));
0484           index++;
0485         }
0486         em_bins.push_back(pt_bin_calibs);
0487       }
0488       l1eg_bins.push_back(em_bins);
0489     }
0490     tauPtCalibrationsHGCal.push_back(l1eg_bins);
0491   }
0492   if (debug) {
0493     LogDebug("L1CaloJetProducer") << " Loading HGCal calibrations: Loaded " << index
0494                                   << " values vs. size() of input calibration file: "
0495                                   << int(tauCalibrationsHGCal.size()) << "\n";
0496   }
0497 
0498   isoTauBarrel.SetParameter(0, 0.30);
0499   isoTauBarrel.SetParameter(1, 0.31);
0500   isoTauBarrel.SetParameter(2, 0.040);
0501   isoTauHGCal.SetParameter(0, 0.34);
0502   isoTauHGCal.SetParameter(1, 0.35);
0503   isoTauHGCal.SetParameter(2, 0.051);
0504 }
0505 
0506 void L1CaloJetProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0507   // Output collections
0508   std::unique_ptr<l1tp2::CaloJetsCollection> L1CaloJetsNoCuts(new l1tp2::CaloJetsCollection);
0509   //std::unique_ptr<l1tp2::CaloJetsCollection> L1CaloJetsWithCuts( new l1tp2::CaloJetsCollection );
0510   //std::unique_ptr<l1extra::L1JetParticleCollection> L1CaloClusterCollectionWithCuts( new l1extra::L1JetParticleCollection );
0511   std::unique_ptr<BXVector<l1t::Jet>> L1CaloJetCollectionBXV(new l1t::JetBxCollection);
0512   std::unique_ptr<BXVector<l1t::Tau>> L1CaloTauCollectionBXV(new l1t::TauBxCollection);
0513 
0514   // Load the ECAL+HCAL tower sums coming from L1EGammaCrystalsEmulatorProducer.cc
0515   std::vector<SimpleCaloHit> l1CaloTowers;
0516 
0517   iEvent.getByToken(l1TowerToken_, l1CaloTowerHandle);
0518   for (auto &hit : *l1CaloTowerHandle.product()) {
0519     SimpleCaloHit l1Hit;
0520     l1Hit.ecalTowerEt = hit.ecalTowerEt();
0521     l1Hit.hcalTowerEt = hit.hcalTowerEt();
0522     l1Hit.l1egTowerEt = hit.l1egTowerEt();
0523     // Add min ET thresholds for tower ET
0524     if (l1Hit.ecalTowerEt < EcalTpEtMin)
0525       l1Hit.ecalTowerEt = 0.0;
0526     if (l1Hit.hcalTowerEt < HcalTpEtMin)
0527       l1Hit.hcalTowerEt = 0.0;
0528     l1Hit.total_tower_et = l1Hit.ecalTowerEt + l1Hit.hcalTowerEt + l1Hit.l1egTowerEt;
0529     l1Hit.towerIEta = hit.towerIEta();
0530     l1Hit.towerIPhi = hit.towerIPhi();
0531     l1Hit.nL1eg = hit.nL1eg();
0532     l1Hit.l1egTrkSS = hit.l1egTrkSS();
0533     l1Hit.l1egTrkIso = hit.l1egTrkIso();
0534     l1Hit.l1egStandaloneSS = hit.l1egStandaloneSS();
0535     l1Hit.l1egStandaloneIso = hit.l1egStandaloneIso();
0536     l1Hit.isBarrel = hit.isBarrel();
0537 
0538     // FIXME There is an error in the L1EGammaCrystalsEmulatorProducer.cc which is
0539     // returning towers with minimal ECAL energy, and no HCAL energy with these
0540     // iEta/iPhi coordinates and eta = -88.653152 and phi = -99.000000.
0541     // Skip these for the time being until the upstream code has been debugged
0542     if ((int)l1Hit.towerIEta == -1016 && (int)l1Hit.towerIPhi == -962)
0543       continue;
0544 
0545     l1Hit.towerEta = hit.towerEta();
0546     l1Hit.towerPhi = hit.towerPhi();
0547     l1CaloTowers.push_back(l1Hit);
0548     if (debug) {
0549       LogDebug("L1CaloJetProducer") << " Tower iEta " << (int)l1Hit.towerIEta << " iPhi " << (int)l1Hit.towerIPhi
0550                                     << " eta " << l1Hit.towerEta << " phi " << l1Hit.towerPhi << " ecal_et "
0551                                     << l1Hit.ecalTowerEt << " hacl_et " << l1Hit.hcalTowerEt << " total_et "
0552                                     << l1Hit.total_tower_et << "\n";
0553     }
0554   }
0555 
0556   // Sort the ECAL+HCAL+L1EGs tower sums based on total ET
0557   std::sort(begin(l1CaloTowers), end(l1CaloTowers), [](const SimpleCaloHit &a, SimpleCaloHit &b) {
0558     return a.total_tower_et > b.total_tower_et;
0559   });
0560 
0561   /**************************************************************************
0562  * Begin with making CaloJets in 9x9 grid based on all energy not included in L1EG Objs.
0563  * For reference, Run-I used 12x12 grid and Stage-2/Phase-I used 9x9 grid.
0564  * We plan to further study this choice and possibly move towards a more circular shape
0565  * Create jetCluster within 9x9 of highest ET seed tower.
0566  * 9 trigger towers contains all of an ak-0.4 jets, but overshoots on the corners.
0567  ******************************************************************************/
0568 
0569   // Experimental parameters, don't want to bother with hardcoding them in data format
0570   std::map<std::string, float> params;
0571 
0572   std::vector<l1CaloJetObj> l1CaloJetObjs;
0573 
0574   // Count the number of unused HCAL TPs so we can stop while loop after done.
0575   // Clustering can also stop once there are no seed hits >= EtMinForSeedHit
0576   int n_towers = l1CaloTowers.size();
0577   int n_stale = 0;
0578   bool caloJetClusteringFinished = false;
0579   while (!caloJetClusteringFinished && n_towers != n_stale) {
0580     l1CaloJetObj caloJetObj;
0581     caloJetObj.Init();
0582 
0583     // First find highest ET ECAL+HCAL+L1EGs tower and use to seed the 9x9 Jet
0584     int cnt = 0;
0585     for (auto &l1CaloTower : l1CaloTowers) {
0586       cnt++;
0587       if (l1CaloTower.stale)
0588         continue;  // skip l1CaloTowers which are already used
0589 
0590       if (caloJetObj.jetClusterET == 0.0)  // this is the first l1CaloTower to seed the jet
0591       {
0592         // Check if the leading unused tower has ET < min for seeding a jet.
0593         // If so, stop jet clustering
0594         if (l1CaloTower.total_tower_et < EtMinForSeedHit) {
0595           caloJetClusteringFinished = true;
0596           continue;
0597         }
0598         l1CaloTower.stale = true;
0599         n_stale++;
0600 
0601         // Set seed location needed for delta iEta/iPhi, eta/phi comparisons later
0602         if (l1CaloTower.isBarrel)
0603           caloJetObj.barrelSeeded = true;
0604         else
0605           caloJetObj.barrelSeeded = false;
0606 
0607         // 3 4-vectors for ECAL, HCAL, ECAL+HCAL for adding together
0608         reco::Candidate::PolarLorentzVector hcalP4(
0609             l1CaloTower.hcalTowerEt, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.);
0610         reco::Candidate::PolarLorentzVector ecalP4(
0611             l1CaloTower.ecalTowerEt, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.);
0612         reco::Candidate::PolarLorentzVector l1egP4(
0613             l1CaloTower.l1egTowerEt, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.);
0614         reco::Candidate::PolarLorentzVector totalP4(
0615             l1CaloTower.total_tower_et, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.);
0616 
0617         if (hcalP4.pt() > 0) {
0618           caloJetObj.hcal_nHits++;
0619           caloJetObj.hcalJetCluster += hcalP4;
0620           caloJetObj.hcalJetClusterET += l1CaloTower.hcalTowerEt;
0621         }
0622         if (ecalP4.pt() > 0) {
0623           caloJetObj.ecal_nHits++;
0624           caloJetObj.ecalJetCluster += ecalP4;
0625           caloJetObj.ecalJetClusterET += l1CaloTower.ecalTowerEt;
0626         }
0627         if (l1egP4.pt() > 0) {
0628           caloJetObj.l1eg_nHits++;
0629           caloJetObj.l1egJetCluster += l1egP4;
0630           caloJetObj.l1egJetClusterET += l1CaloTower.l1egTowerEt;
0631           caloJetObj.l1eg_nL1EGs += l1CaloTower.nL1eg;
0632 
0633           caloJetObj.l1eg_nL1EGs_standaloneSS += l1CaloTower.l1egStandaloneSS;
0634           caloJetObj.l1eg_nL1EGs_standaloneIso += l1CaloTower.l1egStandaloneIso;
0635           caloJetObj.l1eg_nL1EGs_trkMatchSS += l1CaloTower.l1egTrkSS;
0636           caloJetObj.l1eg_nL1EGs_trkMatchIso += l1CaloTower.l1egTrkIso;
0637 
0638           if (l1CaloTower.isBarrel) {
0639             // For decay mode related checks with CaloTaus
0640             // only applicable in the barrel at the moment:
0641             // l1eg pt, HCAL ET, ECAL ET, dEta, dPhi, trkSS, trkIso, standaloneSS, standaloneIso
0642             std::vector<float> l1EG_info = {float(l1egP4.pt()),
0643                                             float(hcalP4.pt()),
0644                                             float(ecalP4.pt()),
0645                                             0.,
0646                                             0.,
0647                                             float(l1CaloTower.l1egTrkSS),
0648                                             float(l1CaloTower.l1egTrkIso),
0649                                             float(l1CaloTower.l1egStandaloneSS),
0650                                             float(l1CaloTower.l1egStandaloneIso)};
0651             if (l1EG_info[1] / (l1EG_info[0] + l1EG_info[2]) < 0.25) {
0652               caloJetObj.n_l1eg_HoverE_LessThreshold++;
0653             }
0654             caloJetObj.associated_l1EGs_.push_back(l1EG_info);
0655           }
0656         }
0657         if (totalP4.pt() > 0) {
0658           caloJetObj.total_nHits++;
0659           caloJetObj.jetCluster += totalP4;
0660           caloJetObj.jetClusterET += l1CaloTower.total_tower_et;
0661           caloJetObj.seedTower += totalP4;
0662           caloJetObj.seedTowerET += l1CaloTower.total_tower_et;
0663         }
0664 
0665         caloJetObj.seed_iEta = l1CaloTower.towerIEta;
0666         caloJetObj.seed_iPhi = l1CaloTower.towerIPhi;
0667 
0668         if (debug) {
0669           LogDebug("L1CaloJetProducer") << " -- hit " << cnt << " , seeding input     p4 pt "
0670                                         << l1CaloTower.total_tower_et << " eta " << l1CaloTower.towerEta << " phi "
0671                                         << l1CaloTower.towerPhi << "\n";
0672           LogDebug("L1CaloJetProducer") << " -- hit " << cnt << " , seeding input2    p4 pt " << totalP4.pt() << " eta "
0673                                         << totalP4.eta() << " phi " << totalP4.phi() << "\n";
0674           LogDebug("L1CaloJetProducer") << " -- hit " << cnt << " seeding reslt tot p4 pt " << caloJetObj.jetClusterET
0675                                         << " eta " << caloJetObj.jetCluster.eta() << " phi "
0676                                         << caloJetObj.jetCluster.phi() << "\n";
0677         }
0678 
0679         // Need to add the seed energy to the dR rings
0680         caloJetObj.hcal_seed += hcalP4.pt();
0681         //caloJetObj.hcal_3x3 += hcalP4.pt();
0682         caloJetObj.hcal_3x5 += hcalP4.pt();
0683         //caloJetObj.hcal_5x5 += hcalP4.pt();
0684         //caloJetObj.hcal_5x7 += hcalP4.pt();
0685         caloJetObj.hcal_7x7 += hcalP4.pt();
0686         caloJetObj.ecal_seed += ecalP4.pt();
0687         //caloJetObj.ecal_3x3 += ecalP4.pt();
0688         caloJetObj.ecal_3x5 += ecalP4.pt();
0689         //caloJetObj.ecal_5x5 += ecalP4.pt();
0690         //caloJetObj.ecal_5x7 += ecalP4.pt();
0691         caloJetObj.ecal_7x7 += ecalP4.pt();
0692         caloJetObj.l1eg_seed += l1egP4.pt();
0693         //caloJetObj.l1eg_3x3 += l1egP4.pt();
0694         caloJetObj.l1eg_3x5 += l1egP4.pt();
0695         //caloJetObj.l1eg_5x5 += l1egP4.pt();
0696         //caloJetObj.l1eg_5x7 += l1egP4.pt();
0697         caloJetObj.l1eg_7x7 += l1egP4.pt();
0698         caloJetObj.total_seed += totalP4.pt();
0699         //caloJetObj.total_3x3 += totalP4.pt();
0700         caloJetObj.total_3x5 += totalP4.pt();
0701         //caloJetObj.total_5x5 += totalP4.pt();
0702         //caloJetObj.total_5x7 += totalP4.pt();
0703         caloJetObj.total_7x7 += totalP4.pt();
0704 
0705         // Some discrimination vars, 2x2s and 2x3 including central seed
0706         //caloJetObj.hcal_2x3 += hcalP4.pt();
0707         //caloJetObj.hcal_2x3_1 += hcalP4.pt();
0708         //caloJetObj.hcal_2x3_2 += hcalP4.pt();
0709         //caloJetObj.ecal_2x3 += ecalP4.pt();
0710         //caloJetObj.ecal_2x3_1 += ecalP4.pt();
0711         //caloJetObj.ecal_2x3_2 += ecalP4.pt();
0712         //caloJetObj.l1eg_2x3 += l1egP4.pt();
0713         //caloJetObj.l1eg_2x3_1 += l1egP4.pt();
0714         //caloJetObj.l1eg_2x3_2 += l1egP4.pt();
0715         //caloJetObj.total_2x3 += totalP4.pt();
0716         //caloJetObj.total_2x3_1 += totalP4.pt();
0717         //caloJetObj.total_2x3_2 += totalP4.pt();
0718 
0719         //caloJetObj.hcal_2x2 += hcalP4.pt();
0720         //caloJetObj.hcal_2x2_1 += hcalP4.pt();
0721         //caloJetObj.hcal_2x2_2 += hcalP4.pt();
0722         //caloJetObj.hcal_2x2_3 += hcalP4.pt();
0723         //caloJetObj.hcal_2x2_4 += hcalP4.pt();
0724         //caloJetObj.ecal_2x2 += ecalP4.pt();
0725         //caloJetObj.ecal_2x2_1 += ecalP4.pt();
0726         //caloJetObj.ecal_2x2_2 += ecalP4.pt();
0727         //caloJetObj.ecal_2x2_3 += ecalP4.pt();
0728         //caloJetObj.ecal_2x2_4 += ecalP4.pt();
0729         //caloJetObj.l1eg_2x2 += l1egP4.pt();
0730         //caloJetObj.l1eg_2x2_1 += l1egP4.pt();
0731         //caloJetObj.l1eg_2x2_2 += l1egP4.pt();
0732         //caloJetObj.l1eg_2x2_3 += l1egP4.pt();
0733         //caloJetObj.l1eg_2x2_4 += l1egP4.pt();
0734         //caloJetObj.total_2x2 += totalP4.pt();
0735         //caloJetObj.total_2x2_1 += totalP4.pt();
0736         //caloJetObj.total_2x2_2 += totalP4.pt();
0737         //caloJetObj.total_2x2_3 += totalP4.pt();
0738         //caloJetObj.total_2x2_4 += totalP4.pt();
0739         continue;
0740       }
0741 
0742       // Unused l1CaloTowers which are not the initial seed
0743       // Depending on seed and tower locations calculate iEta/iPhi or eta/phi comparisons.
0744       // The defaults of 99 will automatically fail comparisons for the incorrect regions.
0745       int hit_iPhi = 99;
0746       int d_iEta = 99;
0747       int d_iPhi = 99;
0748       float d_eta = 99;
0749       float d_phi = 99;
0750       if (caloJetObj.barrelSeeded && l1CaloTower.isBarrel)  // use iEta/iPhi comparisons
0751       {
0752         hit_iPhi = l1CaloTower.towerIPhi;
0753         d_iEta = tower_diEta(caloJetObj.seed_iEta, l1CaloTower.towerIEta);
0754         d_iPhi = tower_diPhi(caloJetObj.seed_iPhi, hit_iPhi);
0755       } else  // either seed or tower are in HGCal or HF, use eta/phi
0756       {
0757         d_eta = caloJetObj.seedTower.eta() - l1CaloTower.towerEta;
0758         d_phi = reco::deltaPhi(caloJetObj.seedTower.phi(), l1CaloTower.towerPhi);
0759       }
0760 
0761       // 7x7 HCAL Trigger Towers
0762       // If seeded in barrel and hit is barrel then we can compare iEta/iPhi, else need to use eta/phi
0763       // in HGCal / transition region
0764       if ((abs(d_iEta) <= 3 && abs(d_iPhi) <= 3) || (std::abs(d_eta) < 0.3 && std::abs(d_phi) < 0.3)) {
0765         l1CaloTower.stale = true;
0766         n_stale++;
0767 
0768         // 3 4-vectors for ECAL, HCAL, ECAL+HCAL for adding together
0769         reco::Candidate::PolarLorentzVector hcalP4(
0770             l1CaloTower.hcalTowerEt, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.);
0771         reco::Candidate::PolarLorentzVector ecalP4(
0772             l1CaloTower.ecalTowerEt, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.);
0773         reco::Candidate::PolarLorentzVector l1egP4(
0774             l1CaloTower.l1egTowerEt, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.);
0775         reco::Candidate::PolarLorentzVector totalP4(
0776             l1CaloTower.total_tower_et, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.);
0777 
0778         if (hcalP4.pt() > 0) {
0779           caloJetObj.hcal_nHits++;
0780           caloJetObj.hcalJetCluster += hcalP4;
0781           caloJetObj.hcalJetClusterET += l1CaloTower.hcalTowerEt;
0782         }
0783         if (ecalP4.pt() > 0) {
0784           caloJetObj.ecal_nHits++;
0785           caloJetObj.ecalJetCluster += ecalP4;
0786           caloJetObj.ecalJetClusterET += l1CaloTower.ecalTowerEt;
0787         }
0788         if (l1egP4.pt() > 0) {
0789           caloJetObj.l1eg_nHits++;
0790           caloJetObj.l1egJetCluster += l1egP4;
0791           caloJetObj.l1egJetClusterET += l1CaloTower.l1egTowerEt;
0792           caloJetObj.l1eg_nL1EGs += l1CaloTower.nL1eg;
0793         }
0794         if (totalP4.pt() > 0) {
0795           caloJetObj.total_nHits++;
0796           caloJetObj.jetCluster += totalP4;
0797           caloJetObj.jetClusterET += l1CaloTower.total_tower_et;
0798         }
0799 
0800         if (debug) {
0801           LogDebug("L1CaloJetProducer") << " ---- hit " << cnt << " input     p4 pt " << totalP4.pt() << " eta "
0802                                         << totalP4.eta() << " phi " << totalP4.phi() << "\n";
0803           LogDebug("L1CaloJetProducer") << " ---- hit " << cnt << " resulting p4 pt " << caloJetObj.jetClusterET
0804                                         << " eta " << caloJetObj.jetCluster.eta() << " phi "
0805                                         << caloJetObj.jetCluster.phi() << "\n";
0806         }
0807 
0808         if ((abs(d_iEta) == 0 && abs(d_iPhi) == 0) || (std::abs(d_eta) < 0.043 && std::abs(d_phi) < 0.043)) {
0809           caloJetObj.hcal_seed += hcalP4.pt();
0810           caloJetObj.ecal_seed += ecalP4.pt();
0811           caloJetObj.l1eg_seed += l1egP4.pt();
0812           caloJetObj.total_seed += totalP4.pt();
0813         }
0814         //if ( (abs( d_iEta ) <= 1 && abs( d_iPhi ) <= 1) ||
0815         //    ( std::abs( d_eta ) < 0.13 && std::abs( d_phi ) < 0.13 ) )
0816         //{
0817         //    caloJetObj.hcal_3x3 += hcalP4.pt();
0818         //    caloJetObj.ecal_3x3 += ecalP4.pt();
0819         //    caloJetObj.l1eg_3x3 += l1egP4.pt();
0820         //    caloJetObj.total_3x3 += totalP4.pt();
0821         //}
0822         if ((abs(d_iEta) <= 1 && abs(d_iPhi) <= 2) || (std::abs(d_eta) < 0.13 && std::abs(d_phi) < 0.22)) {
0823           caloJetObj.hcal_3x5 += hcalP4.pt();
0824           caloJetObj.ecal_3x5 += ecalP4.pt();
0825           caloJetObj.l1eg_3x5 += l1egP4.pt();
0826           caloJetObj.total_3x5 += totalP4.pt();
0827 
0828           // Do this for 3x5 only
0829           if (l1egP4.pt() > 0) {
0830             caloJetObj.l1eg_nL1EGs_standaloneSS += l1CaloTower.l1egStandaloneSS;
0831             caloJetObj.l1eg_nL1EGs_standaloneIso += l1CaloTower.l1egStandaloneIso;
0832             caloJetObj.l1eg_nL1EGs_trkMatchSS += l1CaloTower.l1egTrkSS;
0833             caloJetObj.l1eg_nL1EGs_trkMatchIso += l1CaloTower.l1egTrkIso;
0834 
0835             // For decay mode related checks with CaloTaus
0836             // only applicable in the barrel at the moment:
0837             // l1eg pt, HCAL ET, ECAL ET, d_iEta, d_iPhi, trkSS, trkIso, standaloneSS, standaloneIso
0838             std::vector<float> l1EG_info = {float(l1egP4.pt()),
0839                                             float(hcalP4.pt()),
0840                                             float(ecalP4.pt()),
0841                                             float(d_iEta),
0842                                             float(d_iPhi),
0843                                             float(l1CaloTower.l1egTrkSS),
0844                                             float(l1CaloTower.l1egTrkIso),
0845                                             float(l1CaloTower.l1egStandaloneSS),
0846                                             float(l1CaloTower.l1egStandaloneIso)};
0847             if (l1EG_info[1] / (l1EG_info[0] + l1EG_info[2]) < 0.25) {
0848               caloJetObj.n_l1eg_HoverE_LessThreshold++;
0849             }
0850             caloJetObj.associated_l1EGs_.push_back(l1EG_info);
0851           }
0852         }
0853         //if ( ( abs( d_iEta ) <= 2 && abs( d_iPhi ) <= 2) ||
0854         //    ( std::abs( d_eta ) < 0.22 && std::abs( d_phi ) < 0.22 ) )
0855         //{
0856         //    caloJetObj.hcal_5x5 += hcalP4.pt();
0857         //    caloJetObj.ecal_5x5 += ecalP4.pt();
0858         //    caloJetObj.l1eg_5x5 += l1egP4.pt();
0859         //    caloJetObj.total_5x5 += totalP4.pt();
0860         //}
0861         //if ( ( abs( d_iEta ) <= 2 && abs( d_iPhi ) <= 3) ||
0862         //    ( std::abs( d_eta ) < 0.22 && std::abs( d_phi ) < 0.3 ) )
0863         //{
0864         //    caloJetObj.hcal_5x7 += hcalP4.pt();
0865         //    caloJetObj.ecal_5x7 += ecalP4.pt();
0866         //    caloJetObj.l1eg_5x7 += l1egP4.pt();
0867         //    caloJetObj.total_5x7 += totalP4.pt();
0868         //}
0869         if ((abs(d_iEta) <= 3 && abs(d_iPhi) <= 3) || (std::abs(d_eta) < 0.3 && std::abs(d_phi) < 0.3)) {
0870           caloJetObj.hcal_7x7 += hcalP4.pt();
0871           caloJetObj.ecal_7x7 += ecalP4.pt();
0872           caloJetObj.l1eg_7x7 += l1egP4.pt();
0873           caloJetObj.total_7x7 += totalP4.pt();
0874         }
0875 
0876         //// Some discrimination vars, 2x2s and 2x3s including central seed
0877         //// Barrel first, 2x3
0878         //if ( ( d_iEta == 0 || d_iEta == 1 ) && abs(d_iPhi) <= 1 )
0879         //{
0880         //    caloJetObj.hcal_2x3_1 += hcalP4.pt();
0881         //    caloJetObj.ecal_2x3_1 += ecalP4.pt();
0882         //    caloJetObj.l1eg_2x3_1 += l1egP4.pt();
0883         //    caloJetObj.total_2x3_1 += totalP4.pt();
0884         //}
0885         //if ( ( d_iEta == 0 || d_iEta == -1 ) && abs(d_iPhi) <= 1 )
0886         //{
0887         //    caloJetObj.hcal_2x3_2 += hcalP4.pt();
0888         //    caloJetObj.ecal_2x3_2 += ecalP4.pt();
0889         //    caloJetObj.l1eg_2x3_2 += l1egP4.pt();
0890         //    caloJetObj.total_2x3_2 += totalP4.pt();
0891         //}
0892         //// HGCal / HF
0893         //if ( std::abs( d_eta ) < 0.087 && std::abs( d_phi ) < 0.13 )
0894         //{
0895         //    caloJetObj.hcal_2x3 += hcalP4.pt();
0896         //    caloJetObj.ecal_2x3 += ecalP4.pt();
0897         //    caloJetObj.l1eg_2x3 += l1egP4.pt();
0898         //    caloJetObj.total_2x3 += totalP4.pt();
0899         //}
0900 
0901         //// Now 2x2 Barrel first
0902         //if ( ( d_iEta == 0 || d_iEta == 1 ) && ( d_iPhi == 0 || d_iPhi == 1 ) )
0903         //{
0904         //    caloJetObj.hcal_2x2_1 += hcalP4.pt();
0905         //    caloJetObj.ecal_2x2_1 += ecalP4.pt();
0906         //    caloJetObj.l1eg_2x2_1 += l1egP4.pt();
0907         //    caloJetObj.total_2x2_1 += totalP4.pt();
0908         //}
0909         //if ( ( d_iEta == 0 || d_iEta == 1 ) && ( d_iPhi == 0 || d_iPhi == -1 ) )
0910         //{
0911         //    caloJetObj.hcal_2x2_2 += hcalP4.pt();
0912         //    caloJetObj.ecal_2x2_2 += ecalP4.pt();
0913         //    caloJetObj.l1eg_2x2_2 += l1egP4.pt();
0914         //    caloJetObj.total_2x2_2 += totalP4.pt();
0915         //}
0916         //if ( ( d_iEta == 0 || d_iEta == -1 ) && ( d_iPhi == 0 || d_iPhi == 1 ) )
0917         //{
0918         //    caloJetObj.hcal_2x2_3 += hcalP4.pt();
0919         //    caloJetObj.ecal_2x2_3 += ecalP4.pt();
0920         //    caloJetObj.l1eg_2x2_3 += l1egP4.pt();
0921         //    caloJetObj.total_2x2_3 += totalP4.pt();
0922         //}
0923         //if ( ( d_iEta == 0 || d_iEta == -1 ) && ( d_iPhi == 0 || d_iPhi == -1 ) )
0924         //{
0925         //    caloJetObj.hcal_2x2_4 += hcalP4.pt();
0926         //    caloJetObj.ecal_2x2_4 += ecalP4.pt();
0927         //    caloJetObj.l1eg_2x2_4 += l1egP4.pt();
0928         //    caloJetObj.total_2x2_4 += totalP4.pt();
0929         //}
0930         //// HGCal / HF
0931         //if ( std::abs( d_eta ) < 0.087 && std::abs( d_phi ) < 0.087 )
0932         //{
0933         //    caloJetObj.hcal_2x2 += hcalP4.pt();
0934         //    caloJetObj.ecal_2x2 += ecalP4.pt();
0935         //    caloJetObj.l1eg_2x2 += l1egP4.pt();
0936         //    caloJetObj.total_2x2 += totalP4.pt();
0937         //}
0938       }
0939     }
0940 
0941     if (caloJetObj.jetClusterET > 0.0) {
0942       l1CaloJetObjs.push_back(caloJetObj);
0943     }
0944 
0945   }  // end while loop of HCAL TP clustering
0946 
0947   // Sort JetClusters so we can begin with the highest pt for next step of jet clustering
0948   std::sort(begin(l1CaloJetObjs), end(l1CaloJetObjs), [](const l1CaloJetObj &a, const l1CaloJetObj &b) {
0949     return a.jetClusterET > b.jetClusterET;
0950   });
0951 
0952   /**************************************************************************
0953  * Progress to adding L1EGs built from ECAL TPs  9x9 grid.
0954  * Recall, for 9x9 trigger towers gives diameter 0.78
0955  ******************************************************************************/
0956 
0957   // Cluster together the L1EGs around existing HCAL Jet
0958   // Cluster within dEta/dPhi 0.4 which is very close to 0.39 = 9x9/2
0959   //std::cout << " - Input L1EGs: " << crystalClustersVect.size() << std::endl;
0960   for (auto &caloJetObj : l1CaloJetObjs) {
0961     params["seed_pt"] = caloJetObj.seedTowerET;
0962     params["seed_eta"] = caloJetObj.seedTower.eta();
0963     params["seed_phi"] = caloJetObj.seedTower.phi();
0964     params["seed_iEta"] = caloJetObj.seed_iEta;
0965     params["seed_iPhi"] = caloJetObj.seed_iPhi;
0966     params["seed_energy"] = caloJetObj.seedTower.energy();
0967 
0968     params["hcal_pt"] = caloJetObj.hcalJetClusterET;
0969     params["hcal_seed"] = caloJetObj.hcal_seed;
0970     //params["hcal_3x3"] = caloJetObj.hcal_3x3;
0971     params["hcal_3x5"] = caloJetObj.hcal_3x5;
0972     //params["hcal_5x5"] = caloJetObj.hcal_5x5;
0973     //params["hcal_5x7"] = caloJetObj.hcal_5x7;
0974     params["hcal_7x7"] = caloJetObj.hcal_7x7;
0975     //params["hcal_2x3"] = std::max( caloJetObj.hcal_2x3, std::max( caloJetObj.hcal_2x3_1, caloJetObj.hcal_2x3_2 ));
0976     //params["hcal_2x2"] = std::max( caloJetObj.hcal_2x2, std::max( caloJetObj.hcal_2x2_1, std::max( caloJetObj.hcal_2x2_2, std::max( caloJetObj.hcal_2x2_3, caloJetObj.hcal_2x2_4 ))));
0977     params["hcal_nHits"] = caloJetObj.hcal_nHits;
0978 
0979     params["ecal_pt"] = caloJetObj.ecalJetClusterET;
0980     params["ecal_seed"] = caloJetObj.ecal_seed;
0981     //params["ecal_3x3"] = caloJetObj.ecal_3x3;
0982     params["ecal_3x5"] = caloJetObj.ecal_3x5;
0983     //params["ecal_5x5"] = caloJetObj.ecal_5x5;
0984     //params["ecal_5x7"] = caloJetObj.ecal_5x7;
0985     params["ecal_7x7"] = caloJetObj.ecal_7x7;
0986     //params["ecal_2x3"] = std::max( caloJetObj.ecal_2x3, std::max( caloJetObj.ecal_2x3_1, caloJetObj.ecal_2x3_2 ));
0987     //params["ecal_2x2"] = std::max( caloJetObj.ecal_2x2, std::max( caloJetObj.ecal_2x2_1, std::max( caloJetObj.ecal_2x2_2, std::max( caloJetObj.ecal_2x2_3, caloJetObj.ecal_2x2_4 ))));
0988     params["ecal_nHits"] = caloJetObj.ecal_nHits;
0989 
0990     params["l1eg_pt"] = caloJetObj.l1egJetClusterET;
0991     params["l1eg_seed"] = caloJetObj.l1eg_seed;
0992     //params["l1eg_3x3"] = caloJetObj.l1eg_3x3;
0993     params["l1eg_3x5"] = caloJetObj.l1eg_3x5;
0994     //params["l1eg_5x5"] = caloJetObj.l1eg_5x5;
0995     //params["l1eg_5x7"] = caloJetObj.l1eg_5x7;
0996     params["l1eg_7x7"] = caloJetObj.l1eg_7x7;
0997     //params["l1eg_2x3"] = std::max( caloJetObj.l1eg_2x3, std::max( caloJetObj.l1eg_2x3_1, caloJetObj.l1eg_2x3_2 ));
0998     //params["l1eg_2x2"] = std::max( caloJetObj.l1eg_2x2, std::max( caloJetObj.l1eg_2x2_1, std::max( caloJetObj.l1eg_2x2_2, std::max( caloJetObj.l1eg_2x2_3, caloJetObj.l1eg_2x2_4 ))));
0999     params["l1eg_nHits"] = caloJetObj.l1eg_nHits;
1000     params["l1eg_nL1EGs"] = caloJetObj.l1eg_nL1EGs;
1001     params["l1eg_nL1EGs_standaloneSS"] = caloJetObj.l1eg_nL1EGs_standaloneSS;
1002     params["l1eg_nL1EGs_standaloneIso"] = caloJetObj.l1eg_nL1EGs_standaloneIso;
1003     params["l1eg_nL1EGs_trkMatchSS"] = caloJetObj.l1eg_nL1EGs_trkMatchSS;
1004     params["l1eg_nL1EGs_trkMatchIso"] = caloJetObj.l1eg_nL1EGs_trkMatchIso;
1005 
1006     params["total_et"] = caloJetObj.jetClusterET;
1007     params["total_seed"] = caloJetObj.total_seed;
1008     //params["total_3x3"] = caloJetObj.total_3x3;
1009     params["total_3x5"] = caloJetObj.total_3x5;
1010     //params["total_5x5"] = caloJetObj.total_5x5;
1011     //params["total_5x7"] = caloJetObj.total_5x7;
1012     params["total_7x7"] = caloJetObj.total_7x7;
1013     //params["total_2x3"] = std::max( caloJetObj.total_2x3, std::max( caloJetObj.total_2x3_1, caloJetObj.total_2x3_2 ));
1014     //params["total_2x2"] = std::max( caloJetObj.total_2x2, std::max( caloJetObj.total_2x2_1, std::max( caloJetObj.total_2x2_2, std::max( caloJetObj.total_2x2_3, caloJetObj.total_2x2_4 ))));
1015     params["total_nHits"] = caloJetObj.total_nHits;
1016     //params["total_nTowers"] = total_nTowers;
1017 
1018     //// return -9 for energy and dR values for ecalJet as defaults
1019     float hovere = -9;
1020     if (caloJetObj.ecalJetClusterET > 0.0) {
1021       hovere = caloJetObj.hcalJetClusterET / (caloJetObj.ecalJetClusterET + caloJetObj.l1egJetClusterET);
1022     }
1023 
1024     params["jet_pt"] = caloJetObj.jetClusterET;
1025     params["jet_eta"] = caloJetObj.jetCluster.eta();
1026     params["jet_phi"] = caloJetObj.jetCluster.phi();
1027     params["jet_mass"] = caloJetObj.jetCluster.mass();
1028     params["jet_energy"] = caloJetObj.jetCluster.energy();
1029 
1030     // Calibrations
1031     params["hcal_calibration"] =
1032         get_hcal_calibration(params["jet_pt"], params["ecal_pt"], params["l1eg_pt"], params["jet_eta"]);
1033     params["hcal_pt_calibration"] = params["hcal_pt"] * params["hcal_calibration"];
1034     params["jet_pt_calibration"] = params["hcal_pt_calibration"] + params["ecal_pt"] + params["l1eg_pt"];
1035 
1036     // Tau Vars
1037     // The tau pT calibration is applied as a SF to the total raw pT
1038     // in contrast to the jet calibration above
1039     params["tau_pt_calibration_value"] = get_tau_pt_calibration(params["total_3x5"],
1040                                                                 params["ecal_3x5"],
1041                                                                 params["l1eg_3x5"],
1042                                                                 caloJetObj.n_l1eg_HoverE_LessThreshold,
1043                                                                 params["jet_eta"]);
1044     params["tau_pt"] = params["total_3x5"] * params["tau_pt_calibration_value"];
1045     params["n_l1eg_HoverE_LessThreshold"] = caloJetObj.n_l1eg_HoverE_LessThreshold;
1046     // Currently, applying the tau_pt calibration to the isolation region as well...
1047     // One could switch to using the calibrated jet_pt instead for the iso region...
1048     // This should be revisited - FIXME?
1049     params["tau_total_iso_et"] = params["jet_pt"] * params["tau_pt_calibration_value"];
1050     params["tau_iso_et"] = (params["jet_pt"] * params["tau_pt_calibration_value"]) - params["tau_pt"];
1051     params["loose_iso_tau_wp"] = float(loose_iso_tau_wp(params["tau_pt"], params["tau_iso_et"], params["jet_eta"]));
1052 
1053     float calibratedPt = -1;
1054     float ECalIsolation = -1;  // Need to loop over 7x7 crystals of unclustered energy
1055     float totalPtPUcorr = -1;
1056     l1tp2::CaloJet caloJet(caloJetObj.jetCluster, calibratedPt, hovere, ECalIsolation, totalPtPUcorr);
1057     caloJet.setExperimentalParams(params);
1058     caloJet.setAssociated_l1EGs(caloJetObj.associated_l1EGs_);
1059 
1060     // Only store jets passing ET threshold
1061     if (params["jet_pt_calibration"] >= EtMinForCollection) {
1062       L1CaloJetsNoCuts->push_back(caloJet);
1063       //L1CaloJetsWithCuts->push_back( caloJet );
1064       reco::Candidate::PolarLorentzVector jet_p4 = reco::Candidate::PolarLorentzVector(
1065           params["jet_pt_calibration"], caloJet.p4().eta(), caloJet.p4().phi(), caloJet.p4().M());
1066       L1CaloJetCollectionBXV->push_back(0, l1t::Jet(jet_p4));
1067 
1068       if (debug)
1069         LogDebug("L1CaloJetProducer") << " Made a Jet, eta " << caloJetObj.jetCluster.eta() << " phi "
1070                                       << caloJetObj.jetCluster.phi() << " pt " << caloJetObj.jetClusterET
1071                                       << " calibrated pt " << params["jet_pt_calibration"] << "\n";
1072     }
1073 
1074     // Only store taus passing ET threshold
1075     if (params["tau_pt"] >= EtMinForTauCollection) {
1076       short int tau_ieta = caloJetObj.seed_iEta;
1077       short int tau_iphi = caloJetObj.seed_iPhi;
1078       short int raw_et = params["total_3x5"];
1079       short int iso_et = params["tau_iso_et"];
1080       bool hasEM = false;
1081       if (params["l1eg_3x5"] > 0. || params["ecal_3x5"] > 0.) {
1082         hasEM = true;
1083       }
1084       int tau_qual = int(params["loose_iso_tau_wp"]);
1085 
1086       reco::Candidate::PolarLorentzVector tau_p4 = reco::Candidate::PolarLorentzVector(
1087           params["tau_pt"], caloJet.p4().eta(), caloJet.p4().phi(), caloJet.p4().M());
1088       l1t::Tau l1Tau = l1t::Tau(tau_p4, params["tau_pt"], caloJet.p4().eta(), caloJet.p4().phi(), tau_qual, iso_et);
1089       l1Tau.setTowerIEta(tau_ieta);
1090       l1Tau.setTowerIPhi(tau_iphi);
1091       l1Tau.setRawEt(raw_et);
1092       l1Tau.setIsoEt(iso_et);
1093       l1Tau.setHasEM(hasEM);
1094       l1Tau.setIsMerged(false);
1095       L1CaloTauCollectionBXV->push_back(0, l1Tau);
1096 
1097       if (debug) {
1098         LogDebug("L1CaloJetProducer") << " Made a Jet, eta " << l1Tau.eta() << " phi " << l1Tau.phi() << " pt "
1099                                       << l1Tau.rawEt() << " calibrated pt " << l1Tau.pt() << "\n";
1100       }
1101     }
1102 
1103   }  // end jetClusters loop
1104 
1105   iEvent.put(std::move(L1CaloJetsNoCuts), "L1CaloJetsNoCuts");
1106   //iEvent.put(std::move(L1CaloJetsWithCuts), "L1CaloJetsWithCuts" );
1107   //iEvent.put(std::move(L1CaloClusterCollectionWithCuts), "L1CaloClusterCollectionWithCuts" );
1108   iEvent.put(std::move(L1CaloJetCollectionBXV), "L1CaloJetCollectionBXV");
1109   iEvent.put(std::move(L1CaloTauCollectionBXV), "L1CaloTauCollectionBXV");
1110 }
1111 
1112 int L1CaloJetProducer::tower_diPhi(int &iPhi_1, int &iPhi_2) const {
1113   // 360 Crystals in full, 72 towers, half way is 36
1114   int PI = 36;
1115   int result = iPhi_1 - iPhi_2;
1116   while (result > PI)
1117     result -= 2 * PI;
1118   while (result <= -PI)
1119     result += 2 * PI;
1120   return result;
1121 }
1122 
1123 // Added b/c of the iEta jump from +1 to -1 across the barrel mid point
1124 int L1CaloJetProducer::tower_diEta(int &iEta_1, int &iEta_2) const {
1125   // On same side of barrel
1126   if (iEta_1 * iEta_2 > 0)
1127     return iEta_1 - iEta_2;
1128   else
1129     return iEta_1 - iEta_2 - 1;
1130 }
1131 
1132 float L1CaloJetProducer::get_deltaR(reco::Candidate::PolarLorentzVector &p4_1,
1133                                     reco::Candidate::PolarLorentzVector &p4_2) const {
1134   // Check that pt is > 0 for both or else reco::deltaR returns bogus values
1135   if (p4_1.pt() > 0 && p4_2.pt() > 0) {
1136     return reco::deltaR(p4_1, p4_2);
1137   } else
1138     return -1;
1139 }
1140 
1141 // Apply calibrations to HCAL energy based on EM Fraction, Jet Eta, Jet pT
1142 float L1CaloJetProducer::get_hcal_calibration(float &jet_pt,
1143                                               float &ecal_pt,
1144                                               float &ecal_L1EG_jet_pt,
1145                                               float &jet_eta) const {
1146   float em_frac = (ecal_L1EG_jet_pt + ecal_pt) / jet_pt;
1147   float abs_eta = std::abs(jet_eta);
1148   float tmp_jet_pt = jet_pt;
1149   if (tmp_jet_pt > 499)
1150     tmp_jet_pt = 499;
1151 
1152   // Different indices sizes in different calo regions.
1153   // Barrel...
1154   size_t em_index = 0;
1155   size_t eta_index = 0;
1156   size_t pt_index = 0;
1157   float calib = 1.0;
1158   if (abs_eta <= 1.5) {
1159     // Start loop checking 2nd value
1160     for (unsigned int i = 1; i < emFractionBinsBarrel.size(); i++) {
1161       if (em_frac <= emFractionBinsBarrel.at(i))
1162         break;
1163       em_index++;
1164     }
1165 
1166     // Start loop checking 2nd value
1167     for (unsigned int i = 1; i < absEtaBinsBarrel.size(); i++) {
1168       if (abs_eta <= absEtaBinsBarrel.at(i))
1169         break;
1170       eta_index++;
1171     }
1172 
1173     // Start loop checking 2nd value
1174     for (unsigned int i = 1; i < jetPtBins.size(); i++) {
1175       if (tmp_jet_pt <= jetPtBins.at(i))
1176         break;
1177       pt_index++;
1178     }
1179     calib = calibrationsBarrel[eta_index][em_index][pt_index];
1180   }  // end Barrel
1181   else if (abs_eta <= 3.0)  // HGCal
1182   {
1183     // Start loop checking 2nd value
1184     for (unsigned int i = 1; i < emFractionBinsHGCal.size(); i++) {
1185       if (em_frac <= emFractionBinsHGCal.at(i))
1186         break;
1187       em_index++;
1188     }
1189 
1190     // Start loop checking 2nd value
1191     for (unsigned int i = 1; i < absEtaBinsHGCal.size(); i++) {
1192       if (abs_eta <= absEtaBinsHGCal.at(i))
1193         break;
1194       eta_index++;
1195     }
1196 
1197     // Start loop checking 2nd value
1198     for (unsigned int i = 1; i < jetPtBins.size(); i++) {
1199       if (tmp_jet_pt <= jetPtBins.at(i))
1200         break;
1201       pt_index++;
1202     }
1203     calib = calibrationsHGCal[eta_index][em_index][pt_index];
1204   }  // end HGCal
1205   else  // HF
1206   {
1207     // Start loop checking 2nd value
1208     for (unsigned int i = 1; i < emFractionBinsHF.size(); i++) {
1209       if (em_frac <= emFractionBinsHF.at(i))
1210         break;
1211       em_index++;
1212     }
1213 
1214     // Start loop checking 2nd value
1215     for (unsigned int i = 1; i < absEtaBinsHF.size(); i++) {
1216       if (abs_eta <= absEtaBinsHF.at(i))
1217         break;
1218       eta_index++;
1219     }
1220 
1221     // Start loop checking 2nd value
1222     for (unsigned int i = 1; i < jetPtBins.size(); i++) {
1223       if (tmp_jet_pt <= jetPtBins.at(i))
1224         break;
1225       pt_index++;
1226     }
1227     calib = calibrationsHF[eta_index][em_index][pt_index];
1228   }  // end HF
1229 
1230   return calib;
1231 }
1232 
1233 // Apply calibrations to tau pT based on L1EG info, EM Fraction, Tau Eta, Tau pT
1234 float L1CaloJetProducer::get_tau_pt_calibration(
1235     float &tau_pt, float &ecal_pt, float &l1EG_pt, float &n_L1EGs, float &tau_eta) const {
1236   float em_frac = (l1EG_pt + ecal_pt) / tau_pt;
1237   float abs_eta = std::abs(tau_eta);
1238   float tmp_tau_pt = tau_pt;
1239   if (tmp_tau_pt > 199)
1240     tmp_tau_pt = 199;
1241 
1242   // Different indices sizes in different calo regions.
1243   // Barrel...
1244   size_t em_index = 0;
1245   size_t eta_index = 0;
1246   size_t n_L1EG_index = 0;
1247   size_t pt_index = 0;
1248   float calib = 1.0;
1249   // HERE
1250   if (abs_eta <= 1.5) {
1251     // Start loop checking 1st value
1252     for (unsigned int i = 0; i < tauL1egValuesBarrel.size(); i++) {
1253       if (n_L1EGs == tauL1egValuesBarrel.at(i))
1254         break;
1255       if (tauL1egValuesBarrel.at(i) == tauL1egValuesBarrel.back())
1256         break;  // to preven incrementing on last one
1257       n_L1EG_index++;
1258     }
1259 
1260     // Find key value pair matching n L1EGs
1261     for (auto &l1eg_info : tauL1egInfoMapBarrel) {
1262       if (l1eg_info.first != double(n_L1EG_index))
1263         continue;
1264       // Start loop checking 2nd value
1265       for (unsigned int i = 1; i < l1eg_info.second.size(); i++) {
1266         if (em_frac <= l1eg_info.second.at(i))
1267           break;
1268         em_index++;
1269       }
1270     }
1271 
1272     // Start loop checking 2nd value
1273     for (unsigned int i = 1; i < tauAbsEtaBinsBarrel.size(); i++) {
1274       if (abs_eta <= tauAbsEtaBinsBarrel.at(i))
1275         break;
1276       eta_index++;
1277     }
1278 
1279     // Start loop checking 2nd value
1280     for (unsigned int i = 1; i < tauPtBins.size(); i++) {
1281       if (tmp_tau_pt <= tauPtBins.at(i))
1282         break;
1283       pt_index++;
1284     }
1285     calib = tauPtCalibrationsBarrel[eta_index][n_L1EG_index][em_index][pt_index];
1286   }  // end Barrel
1287   else if (abs_eta <= 3.0)  // HGCal
1288   {
1289     // Start loop checking 1st value
1290     for (unsigned int i = 0; i < tauL1egValuesHGCal.size(); i++) {
1291       if (n_L1EGs == tauL1egValuesHGCal.at(i))
1292         break;
1293       if (tauL1egValuesHGCal.at(i) == tauL1egValuesHGCal.back())
1294         break;  // to preven incrementing on last one
1295       n_L1EG_index++;
1296     }
1297 
1298     // Find key value pair matching n L1EGs
1299     for (auto &l1eg_info : tauL1egInfoMapHGCal) {
1300       if (l1eg_info.first != double(n_L1EG_index))
1301         continue;
1302       // Start loop checking 2nd value
1303       for (unsigned int i = 1; i < l1eg_info.second.size(); i++) {
1304         if (em_frac <= l1eg_info.second.at(i))
1305           break;
1306         em_index++;
1307       }
1308     }
1309 
1310     // Start loop checking 2nd value
1311     for (unsigned int i = 1; i < tauAbsEtaBinsHGCal.size(); i++) {
1312       if (abs_eta <= tauAbsEtaBinsHGCal.at(i))
1313         break;
1314       eta_index++;
1315     }
1316 
1317     // Start loop checking 2nd value
1318     for (unsigned int i = 1; i < tauPtBins.size(); i++) {
1319       if (tmp_tau_pt <= tauPtBins.at(i))
1320         break;
1321       pt_index++;
1322     }
1323     calib = tauPtCalibrationsHGCal[eta_index][n_L1EG_index][em_index][pt_index];
1324   }  // end HGCal
1325   else
1326     return calib;
1327 
1328   return calib;
1329 }
1330 
1331 // Loose IsoTau WP
1332 int L1CaloJetProducer::loose_iso_tau_wp(float &tau_pt, float &tau_iso_et, float &tau_eta) const {
1333   // Fully relaxed above 100 GeV pT
1334   if (tau_pt > 100) {
1335     return 1;
1336   }
1337   // Split by barrel and HGCal
1338   // with Barrel first
1339   if (std::abs(tau_eta) < 1.5) {
1340     if (isoTauBarrel.Eval(tau_pt) >= (tau_iso_et / tau_pt)) {
1341       return 1;
1342     } else {
1343       return 0;
1344     }
1345   }
1346   // HGCal
1347   if (std::abs(tau_eta) < 3.0) {
1348     if (isoTauHGCal.Eval(tau_pt) >= (tau_iso_et / tau_pt)) {
1349       return 1;
1350     } else {
1351       return 0;
1352     }
1353   }
1354   // Beyond HGCal
1355   return 0;
1356 }
1357 
1358 DEFINE_FWK_MODULE(L1CaloJetProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!