Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:20:10 

0001 // -*- C++ -*-
0002 //
0003 // Package: L1CaloTrigger
0004 // Class: L1TowerCalibrator
0005 //
0006 /**\class L1TowerCalibrator L1TowerCalibrator.cc
0007 
0008 Description: 
0009 Take the calibrated unclustered ECAL energy and total HCAL
0010 energy associated with the L1CaloTower collection output
0011 from L1EGammaCrystalsEmulatorProducer: l1CaloTowerCollection, "L1CaloTowerCollection"
0012 
0013 as well as HGCal Tower level inputs:
0014 BXVector<l1t::HGCalTower> "hgcalTriggerPrimitiveDigiProducer" "tower" "HLT"     
0015 
0016 and HCAL HF inputs from:
0017 edm::SortedCollection<HcalTriggerPrimitiveDigi,edm::StrictWeakOrdering<HcalTriggerPrimitiveDigi> > "simHcalTriggerPrimitiveDigis" "" "HLT"     
0018 
0019 
0020 Implement PU-based calibrations which scale down the ET
0021 in the towers based on mapping nTowers with ECAL(HCAL) ET <= defined PU threshold.
0022 This value has been shown to be similar between TTbar, QCD, and minBias samples.
0023 This allows a prediction of nvtx. Which can be mapped to the total minBias
0024 energy in an eta slice of the detector.  Subtract the total estimated minBias energy
0025 per eta slice divided by nTowers in that eta slice from each trigger tower in
0026 that eta slice.
0027 
0028 This is all ECAL / HCAL specific or EM vs. Hadronic for HGCal.
0029 
0030 Implementation:
0031 [Notes on implementation]
0032 */
0033 //
0034 // Original Author: Tyler Ruggles
0035 // Created: Thu Nov 15 2018
0036 // $Id$
0037 //
0038 //
0039 
0040 #include "FWCore/Framework/interface/Frameworkfwd.h"
0041 #include "FWCore/Framework/interface/one/EDProducer.h"
0042 #include "FWCore/Framework/interface/ESHandle.h"
0043 #include "FWCore/Framework/interface/Event.h"
0044 #include "FWCore/Framework/interface/MakerMacros.h"
0045 #include "FWCore/ServiceRegistry/interface/Service.h"
0046 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0047 
0048 #include <iostream>
0049 
0050 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloTower.h"
0051 #include "DataFormats/L1THGCal/interface/HGCalTower.h"
0052 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0053 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0054 
0055 #include "CalibFormats/CaloTPG/interface/CaloTPGTranscoder.h"
0056 #include "CalibFormats/CaloTPG/interface/CaloTPGRecord.h"
0057 
0058 #include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
0059 
0060 #include "TFile.h"
0061 #include "TF1.h"
0062 
0063 class L1TowerCalibrator : public edm::one::EDProducer<> {
0064 public:
0065   explicit L1TowerCalibrator(const edm::ParameterSet&);
0066 
0067 private:
0068   void produce(edm::Event&, const edm::EventSetup&) override;
0069 
0070   const double HcalTpEtMin;
0071   const double EcalTpEtMin;
0072   const double HGCalHadTpEtMin;
0073   const double HGCalEmTpEtMin;
0074   const double HFTpEtMin;
0075   const double puThreshold;
0076   const double puThresholdL1eg;
0077   const double puThresholdHcalMin;
0078   const double puThresholdHcalMax;
0079   const double puThresholdEcalMin;
0080   const double puThresholdEcalMax;
0081   const double puThresholdHGCalEMMin;
0082   const double puThresholdHGCalEMMax;
0083   const double puThresholdHGCalHadMin;
0084   const double puThresholdHGCalHadMax;
0085   const double puThresholdHFMin;
0086   const double puThresholdHFMax;
0087   const double barrelSF;
0088   const double hgcalSF;
0089   const double hfSF;
0090   const bool debug;
0091   const bool skipCalibrations;
0092 
0093   const edm::EDGetTokenT<l1tp2::CaloTowerCollection> l1TowerToken_;
0094   edm::Handle<l1tp2::CaloTowerCollection> l1CaloTowerHandle;
0095 
0096   const edm::EDGetTokenT<l1t::HGCalTowerBxCollection> hgcalTowersToken_;
0097   edm::Handle<l1t::HGCalTowerBxCollection> hgcalTowersHandle;
0098   l1t::HGCalTowerBxCollection hgcalTowers;
0099 
0100   const edm::EDGetTokenT<HcalTrigPrimDigiCollection> hcalToken_;
0101   edm::Handle<HcalTrigPrimDigiCollection> hcalTowerHandle;
0102   const edm::ESGetToken<CaloTPGTranscoder, CaloTPGRecord> decoderTag_;
0103 
0104   // nHits to nvtx functions
0105   std::vector<edm::ParameterSet> nHits_to_nvtx_params;
0106   std::map<std::string, TF1> nHits_to_nvtx_funcs;
0107 
0108   // nvtx to PU subtraction functions
0109   std::vector<edm::ParameterSet> nvtx_to_PU_sub_params;
0110   std::map<std::string, TF1> ecal_nvtx_to_PU_sub_funcs;
0111   std::map<std::string, TF1> hcal_nvtx_to_PU_sub_funcs;
0112   std::map<std::string, TF1> hgcalEM_nvtx_to_PU_sub_funcs;
0113   std::map<std::string, TF1> hgcalHad_nvtx_to_PU_sub_funcs;
0114   std::map<std::string, TF1> hf_nvtx_to_PU_sub_funcs;
0115   std::map<std::string, std::map<std::string, TF1> > all_nvtx_to_PU_sub_funcs;
0116 };
0117 
0118 L1TowerCalibrator::L1TowerCalibrator(const edm::ParameterSet& iConfig)
0119     : HcalTpEtMin(iConfig.getParameter<double>("HcalTpEtMin")),
0120       EcalTpEtMin(iConfig.getParameter<double>("EcalTpEtMin")),
0121       HGCalHadTpEtMin(iConfig.getParameter<double>("HGCalHadTpEtMin")),
0122       HGCalEmTpEtMin(iConfig.getParameter<double>("HGCalEmTpEtMin")),
0123       HFTpEtMin(iConfig.getParameter<double>("HFTpEtMin")),
0124       puThreshold(iConfig.getParameter<double>("puThreshold")),
0125       puThresholdL1eg(iConfig.getParameter<double>("puThresholdL1eg")),
0126       puThresholdHcalMin(iConfig.getParameter<double>("puThresholdHcalMin")),
0127       puThresholdHcalMax(iConfig.getParameter<double>("puThresholdHcalMax")),
0128       puThresholdEcalMin(iConfig.getParameter<double>("puThresholdEcalMin")),
0129       puThresholdEcalMax(iConfig.getParameter<double>("puThresholdEcalMax")),
0130       puThresholdHGCalEMMin(iConfig.getParameter<double>("puThresholdHGCalEMMin")),
0131       puThresholdHGCalEMMax(iConfig.getParameter<double>("puThresholdHGCalEMMax")),
0132       puThresholdHGCalHadMin(iConfig.getParameter<double>("puThresholdHGCalHadMin")),
0133       puThresholdHGCalHadMax(iConfig.getParameter<double>("puThresholdHGCalHadMax")),
0134       puThresholdHFMin(iConfig.getParameter<double>("puThresholdHFMin")),
0135       puThresholdHFMax(iConfig.getParameter<double>("puThresholdHFMax")),
0136       barrelSF(iConfig.getParameter<double>("barrelSF")),
0137       hgcalSF(iConfig.getParameter<double>("hgcalSF")),
0138       hfSF(iConfig.getParameter<double>("hfSF")),
0139       debug(iConfig.getParameter<bool>("debug")),
0140       skipCalibrations(iConfig.getParameter<bool>("skipCalibrations")),
0141       l1TowerToken_(consumes<l1tp2::CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("l1CaloTowers"))),
0142       hgcalTowersToken_(
0143           consumes<l1t::HGCalTowerBxCollection>(iConfig.getParameter<edm::InputTag>("L1HgcalTowersInputTag"))),
0144       hcalToken_(consumes<HcalTrigPrimDigiCollection>(iConfig.getParameter<edm::InputTag>("hcalDigis"))),
0145       decoderTag_(esConsumes<CaloTPGTranscoder, CaloTPGRecord>(edm::ESInputTag("", ""))),
0146       nHits_to_nvtx_params(iConfig.getParameter<std::vector<edm::ParameterSet> >("nHits_to_nvtx_params")),
0147       nvtx_to_PU_sub_params(iConfig.getParameter<std::vector<edm::ParameterSet> >("nvtx_to_PU_sub_params")) {
0148   // Initialize the nHits --> nvtx functions
0149   for (uint i = 0; i < nHits_to_nvtx_params.size(); i++) {
0150     edm::ParameterSet* pset = &nHits_to_nvtx_params.at(i);
0151     std::string calo = pset->getParameter<std::string>("fit");
0152     nHits_to_nvtx_funcs[calo.c_str()] = TF1(calo.c_str(), "[0] + [1] * x");
0153     nHits_to_nvtx_funcs[calo.c_str()].SetParameter(0, pset->getParameter<std::vector<double> >("params").at(0));
0154     nHits_to_nvtx_funcs[calo.c_str()].SetParameter(1, pset->getParameter<std::vector<double> >("params").at(1));
0155 
0156     if (debug) {
0157       printf(
0158           "nHits_to_nvtx_params[%i]\n \
0159                 fit: %s \n \
0160                 p1: %f \n \
0161                 p2 %f \n",
0162           i,
0163           calo.c_str(),
0164           nHits_to_nvtx_funcs[calo.c_str()].GetParameter(0),
0165           nHits_to_nvtx_funcs[calo.c_str()].GetParameter(1));
0166     }
0167   }
0168 
0169   // Initialize the nvtx --> PU subtraction functions
0170   all_nvtx_to_PU_sub_funcs["ecal"] = ecal_nvtx_to_PU_sub_funcs;
0171   all_nvtx_to_PU_sub_funcs["hcal"] = hcal_nvtx_to_PU_sub_funcs;
0172   all_nvtx_to_PU_sub_funcs["hgcalEM"] = hgcalEM_nvtx_to_PU_sub_funcs;
0173   all_nvtx_to_PU_sub_funcs["hgcalHad"] = hgcalHad_nvtx_to_PU_sub_funcs;
0174   all_nvtx_to_PU_sub_funcs["hf"] = hf_nvtx_to_PU_sub_funcs;
0175 
0176   for (uint i = 0; i < nvtx_to_PU_sub_params.size(); i++) {
0177     edm::ParameterSet* pset = &nvtx_to_PU_sub_params.at(i);
0178     std::string calo = pset->getParameter<std::string>("calo");
0179     std::string iEta = pset->getParameter<std::string>("iEta");
0180     double p1 = pset->getParameter<std::vector<double> >("params").at(0);
0181     double p2 = pset->getParameter<std::vector<double> >("params").at(1);
0182 
0183     all_nvtx_to_PU_sub_funcs[calo.c_str()][iEta.c_str()] = TF1(calo.c_str(), "[0] + [1] * x");
0184     all_nvtx_to_PU_sub_funcs[calo.c_str()][iEta.c_str()].SetParameter(0, p1);
0185     all_nvtx_to_PU_sub_funcs[calo.c_str()][iEta.c_str()].SetParameter(1, p2);
0186 
0187     if (debug) {
0188       printf(
0189           "nvtx_to_PU_sub_params[%i]\n \
0190                 sub detector: %s \n \
0191                 iEta: %s \n \
0192                 p1: %f \n \
0193                 p2 %f \n",
0194           i,
0195           calo.c_str(),
0196           iEta.c_str(),
0197           p1,
0198           p2);
0199     }
0200   }
0201 
0202   // Our two outputs, calibrated towers and estimated nvtx for fun
0203   produces<l1tp2::CaloTowerCollection>("L1CaloTowerCalibratedCollection");
0204   produces<double>("EstimatedNvtx");
0205 }
0206 
0207 void L1TowerCalibrator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0208   // Estimated number of vertices used for calibration estimattion
0209   std::unique_ptr<double> EstimatedNvtx(new double);
0210   // Calibrated output collection
0211   std::unique_ptr<l1tp2::CaloTowerCollection> L1CaloTowerCalibratedCollection(new l1tp2::CaloTowerCollection);
0212 
0213   // Load the ECAL+HCAL tower sums coming from L1EGammaCrystalsEmulatorProducer.cc
0214   iEvent.getByToken(l1TowerToken_, l1CaloTowerHandle);
0215 
0216   // HGCal info
0217   iEvent.getByToken(hgcalTowersToken_, hgcalTowersHandle);
0218   hgcalTowers = (*hgcalTowersHandle.product());
0219 
0220   // HF Tower info
0221   iEvent.getByToken(hcalToken_, hcalTowerHandle);
0222 
0223   // Barrel ECAL (unclustered) and HCAL
0224   for (auto& hit : *l1CaloTowerHandle.product()) {
0225     l1tp2::CaloTower l1Hit;
0226     l1Hit.setEcalTowerEt(hit.ecalTowerEt());
0227     l1Hit.setHcalTowerEt(hit.hcalTowerEt());
0228     l1Hit.setL1egTowerEt(hit.l1egTowerEt());
0229     // Add min ET thresholds for tower ET
0230     if (l1Hit.ecalTowerEt() < EcalTpEtMin)
0231       l1Hit.setEcalTowerEt(0.0);
0232     if (l1Hit.hcalTowerEt() < HcalTpEtMin)
0233       l1Hit.setHcalTowerEt(0.0);
0234     l1Hit.setTowerIEta(hit.towerIEta());
0235     l1Hit.setTowerIPhi(hit.towerIPhi());
0236     l1Hit.setTowerEta(hit.towerEta());
0237     l1Hit.setTowerPhi(hit.towerPhi());
0238     l1Hit.setIsBarrel(hit.isBarrel());
0239     l1Hit.setNL1eg(hit.nL1eg());
0240     l1Hit.setL1egTrkSS(hit.l1egTrkSS());
0241     l1Hit.setL1egTrkIso(hit.l1egTrkIso());
0242     l1Hit.setL1egStandaloneSS(hit.l1egStandaloneSS());
0243     l1Hit.setL1egStandaloneIso(hit.l1egStandaloneIso());
0244 
0245     // FIXME There is an error in the L1EGammaCrystalsEmulatorProducer.cc which is
0246     // returning towers with minimal ECAL energy, and no HCAL energy with these
0247     // iEta/iPhi coordinates and eta = -88.653152 and phi = -99.000000.
0248     // Skip these for the time being until the upstream code has been debugged
0249     if ((int)l1Hit.towerIEta() == -1016 && (int)l1Hit.towerIPhi() == -962)
0250       continue;
0251 
0252     (*L1CaloTowerCalibratedCollection).push_back(l1Hit);
0253     if (debug)
0254       printf("Barrel tower iEta %i iPhi %i eta %f phi %f ecal_et %f hcal_et_sum %f\n",
0255              (int)l1Hit.towerIEta(),
0256              (int)l1Hit.towerIPhi(),
0257              l1Hit.towerEta(),
0258              l1Hit.towerPhi(),
0259              l1Hit.ecalTowerEt(),
0260              l1Hit.hcalTowerEt());
0261   }
0262 
0263   // Loop over HGCalTowers and create L1CaloTowers for them and add to collection
0264   // This iterator is taken from the PF P2 group
0265   // https://github.com/p2l1pfp/cmssw/blob/170808db68038d53794bc65fdc962f8fc337a24d/L1Trigger/Phase2L1ParticleFlow/plugins/L1TPFCaloProducer.cc#L278-L289
0266   for (auto it = hgcalTowers.begin(0), ed = hgcalTowers.end(0); it != ed; ++it) {
0267     // skip lowest ET towers
0268     if (it->etEm() < HGCalEmTpEtMin && it->etHad() < HGCalHadTpEtMin)
0269       continue;
0270 
0271     l1tp2::CaloTower l1Hit;
0272     // Set energies normally, but need to zero if below threshold
0273     if (it->etEm() < HGCalEmTpEtMin)
0274       l1Hit.setEcalTowerEt(0.);
0275     else
0276       l1Hit.setEcalTowerEt(it->etEm());
0277 
0278     if (it->etHad() < HGCalHadTpEtMin)
0279       l1Hit.setHcalTowerEt(0.);
0280     else
0281       l1Hit.setHcalTowerEt(it->etHad());
0282 
0283     l1Hit.setTowerEta(it->eta());
0284     l1Hit.setTowerPhi(it->phi());
0285     l1Hit.setTowerIEta(-98);  // -98 mean HGCal
0286     l1Hit.setTowerIPhi(-98);
0287     l1Hit.setIsBarrel(false);
0288     (*L1CaloTowerCalibratedCollection).push_back(l1Hit);
0289     if (debug)
0290       printf("HGCal tower iEta %i iPhi %i eta %f phi %f ecal_et %f hcal_et_sum %f\n",
0291              (int)l1Hit.towerIEta(),
0292              (int)l1Hit.towerIPhi(),
0293              l1Hit.towerEta(),
0294              l1Hit.towerPhi(),
0295              l1Hit.ecalTowerEt(),
0296              l1Hit.hcalTowerEt());
0297   }
0298 
0299   // Loop over Hcal HF tower inputs and create L1CaloTowers and add to
0300   // L1CaloTowerCalibratedCollection collection
0301   const auto& decoder = iSetup.getData(decoderTag_);
0302   for (const auto& hit : *hcalTowerHandle.product()) {
0303     HcalTrigTowerDetId id = hit.id();
0304     double et = decoder.hcaletValue(hit.id(), hit.t0());
0305     if (et < HFTpEtMin)
0306       continue;
0307     // Only doing HF so skip outside range
0308     if (abs(id.ieta()) < l1t::CaloTools::kHFBegin)
0309       continue;
0310     if (abs(id.ieta()) > l1t::CaloTools::kHFEnd)
0311       continue;
0312 
0313     l1tp2::CaloTower l1Hit;
0314     l1Hit.setEcalTowerEt(0.);
0315     l1Hit.setHcalTowerEt(et);
0316     l1Hit.setTowerEta(l1t::CaloTools::towerEta(id.ieta()));
0317     l1Hit.setTowerPhi(l1t::CaloTools::towerPhi(id.ieta(), id.iphi()));
0318     l1Hit.setTowerIEta(id.ieta());
0319     l1Hit.setTowerIPhi(id.iphi());
0320     l1Hit.setIsBarrel(false);
0321     (*L1CaloTowerCalibratedCollection).push_back(l1Hit);
0322 
0323     if (debug)
0324       printf("HCAL HF tower iEta %i iPhi %i eta %f phi %f ecal_et %f hcal_et_sum %f\n",
0325              (int)l1Hit.towerIEta(),
0326              (int)l1Hit.towerIPhi(),
0327              l1Hit.towerEta(),
0328              l1Hit.towerPhi(),
0329              l1Hit.ecalTowerEt(),
0330              l1Hit.hcalTowerEt());
0331   }
0332 
0333   // N Tower totals
0334   // For mapping to estimated nvtx in event
0335   int i_ecal_hits_leq_threshold = 0;
0336   int i_hgcalEM_hits_leq_threshold = 0;
0337   int i_hcal_hits_leq_threshold = 0;
0338   int i_hgcalHad_hits_leq_threshold = 0;
0339   int i_hf_hits_leq_threshold = 0;
0340 
0341   // Loop over the collection containing all hits
0342   // and calculate the number of hits falling into the
0343   // "less than or equal" nTowers variable which maps to
0344   // estimated number of vertices
0345   for (auto& l1CaloTower : (*L1CaloTowerCalibratedCollection)) {
0346     // Barrel ECAL
0347     if (l1CaloTower.ecalTowerEt() > 0. && l1CaloTower.towerIEta() != -98) {
0348       if (l1CaloTower.ecalTowerEt() <= puThresholdEcalMax && l1CaloTower.ecalTowerEt() >= puThresholdEcalMin) {
0349         i_ecal_hits_leq_threshold++;
0350       }
0351     }
0352 
0353     // HGCal EM
0354     if (l1CaloTower.ecalTowerEt() > 0. && l1CaloTower.towerIEta() == -98) {
0355       if (l1CaloTower.ecalTowerEt() <= puThresholdHGCalEMMax && l1CaloTower.ecalTowerEt() >= puThresholdHGCalEMMin) {
0356         i_hgcalEM_hits_leq_threshold++;
0357       }
0358     }
0359 
0360     // Barrel HCAL
0361     if (l1CaloTower.hcalTowerEt() > 0. && l1CaloTower.towerIEta() != -98 &&
0362         abs(l1CaloTower.towerEta()) < 2.0)  // abs(eta) < 2 just keeps us out of HF
0363     {
0364       if (l1CaloTower.hcalTowerEt() <= puThresholdHcalMax && l1CaloTower.hcalTowerEt() >= puThresholdHcalMin) {
0365         i_hcal_hits_leq_threshold++;
0366       }
0367     }
0368 
0369     // HGCal Had
0370     if (l1CaloTower.hcalTowerEt() > 0. && l1CaloTower.towerIEta() == -98) {
0371       if (l1CaloTower.hcalTowerEt() <= puThresholdHGCalHadMax && l1CaloTower.hcalTowerEt() >= puThresholdHGCalHadMin) {
0372         i_hgcalHad_hits_leq_threshold++;
0373       }
0374     }
0375 
0376     // HF
0377     if (l1CaloTower.hcalTowerEt() > 0. && l1CaloTower.towerIEta() != -98 &&
0378         abs(l1CaloTower.towerEta()) > 2.0)  // abs(eta) > 2 keeps us out of barrel HF
0379     {
0380       if (l1CaloTower.hcalTowerEt() <= puThresholdHFMax && l1CaloTower.hcalTowerEt() >= puThresholdHFMin) {
0381         i_hf_hits_leq_threshold++;
0382       }
0383     }
0384   }
0385 
0386   // For each subdetector, map to the estimated number of vertices
0387   double ecal_nvtx = nHits_to_nvtx_funcs["ecal"].Eval(i_ecal_hits_leq_threshold);
0388   double hcal_nvtx = nHits_to_nvtx_funcs["hcal"].Eval(i_hcal_hits_leq_threshold);
0389   double hgcalEM_nvtx = nHits_to_nvtx_funcs["hgcalEM"].Eval(i_hgcalEM_hits_leq_threshold);
0390   double hgcalHad_nvtx = nHits_to_nvtx_funcs["hgcalHad"].Eval(i_hgcalHad_hits_leq_threshold);
0391   double hf_nvtx = nHits_to_nvtx_funcs["hf"].Eval(i_hf_hits_leq_threshold);
0392   // Make sure all values are >= 0
0393   if (ecal_nvtx < 0)
0394     ecal_nvtx = 0;
0395   if (hcal_nvtx < 0)
0396     hcal_nvtx = 0;
0397   if (hgcalEM_nvtx < 0)
0398     hgcalEM_nvtx = 0;
0399   if (hgcalHad_nvtx < 0)
0400     hgcalHad_nvtx = 0;
0401   if (hf_nvtx < 0)
0402     hf_nvtx = 0;
0403   // Best estimate is avg of all except HF.
0404   // This is b/c HF currently has such poor prediction power, it only degrades the avg result
0405   // NEW, with 10_3_X, hgcal and HF has the best results based on the values I took...
0406   // skip ECAL and HCAL
0407   //*EstimatedNvtx = ( ecal_nvtx + hcal_nvtx + hgcalEM_nvtx + hgcalHad_nvtx + hf_nvtx ) / 3.;
0408   *EstimatedNvtx = (hgcalEM_nvtx + hgcalHad_nvtx + hf_nvtx) / 3.;
0409 
0410   if (debug) {
0411     double lumi = iEvent.eventAuxiliary().luminosityBlock();
0412     double event = iEvent.eventAuxiliary().event();
0413 
0414     printf(
0415         "L1TowerCalibrater: lumi %.0f evt %.0f nTowers for subdetecters \
0416             \nECAL:      %i --> nvtx = %.1f \
0417             \nHGCal EM:  %i --> nvtx = %.1f \
0418             \nHCAL:      %i --> nvtx = %.1f \
0419             \nHGCal Had: %i --> nvtx = %.1f \
0420             \nHCAL HF:   %i --> nvtx = %.1f \
0421             \nEstimated Nvtx = %.1f\n",
0422         lumi,
0423         event,
0424         i_ecal_hits_leq_threshold,
0425         ecal_nvtx,
0426         i_hgcalEM_hits_leq_threshold,
0427         hgcalEM_nvtx,
0428         i_hcal_hits_leq_threshold,
0429         hcal_nvtx,
0430         i_hgcalHad_hits_leq_threshold,
0431         hgcalHad_nvtx,
0432         i_hf_hits_leq_threshold,
0433         hf_nvtx,
0434         *EstimatedNvtx);
0435   }
0436 
0437   // Use estimated number of vertices to subtract off PU contributions
0438   // to each and every hit. In cases where the energy would go negative,
0439   // limit this to zero.
0440   if (!skipCalibrations)  // skipCalibrations simply passes the towers through
0441   {
0442     for (auto& l1CaloTower : (*L1CaloTowerCalibratedCollection)) {
0443       // Barrel ECAL eta slices
0444       if (l1CaloTower.ecalTowerEt() > 0. && l1CaloTower.towerIEta() != -98) {
0445         if (abs(l1CaloTower.towerIEta()) <= 3) {
0446           l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() -
0447                                      all_nvtx_to_PU_sub_funcs["ecal"]["er1to3"].Eval(*EstimatedNvtx) * barrelSF);
0448         }
0449         if (abs(l1CaloTower.towerIEta()) <= 6 && abs(l1CaloTower.towerIEta()) >= 4) {
0450           l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() -
0451                                      all_nvtx_to_PU_sub_funcs["ecal"]["er4to6"].Eval(*EstimatedNvtx) * barrelSF);
0452         }
0453         if (abs(l1CaloTower.towerIEta()) <= 9 && abs(l1CaloTower.towerIEta()) >= 7) {
0454           l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() -
0455                                      all_nvtx_to_PU_sub_funcs["ecal"]["er7to9"].Eval(*EstimatedNvtx) * barrelSF);
0456         }
0457         if (abs(l1CaloTower.towerIEta()) <= 12 && abs(l1CaloTower.towerIEta()) >= 10) {
0458           l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() -
0459                                      all_nvtx_to_PU_sub_funcs["ecal"]["er10to12"].Eval(*EstimatedNvtx) * barrelSF);
0460         }
0461         if (abs(l1CaloTower.towerIEta()) <= 15 && abs(l1CaloTower.towerIEta()) >= 13) {
0462           l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() -
0463                                      all_nvtx_to_PU_sub_funcs["ecal"]["er13to15"].Eval(*EstimatedNvtx) * barrelSF);
0464         }
0465         if (abs(l1CaloTower.towerIEta()) <= 18 && abs(l1CaloTower.towerIEta()) >= 16) {
0466           l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() -
0467                                      all_nvtx_to_PU_sub_funcs["ecal"]["er16to18"].Eval(*EstimatedNvtx) * barrelSF);
0468         }
0469       }
0470 
0471       // HGCal EM eta slices
0472       if (l1CaloTower.ecalTowerEt() > 0. && l1CaloTower.towerIEta() == -98) {
0473         if (abs(l1CaloTower.towerEta()) <= 1.8) {
0474           l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() -
0475                                      all_nvtx_to_PU_sub_funcs["hgcalEM"]["er1p4to1p8"].Eval(*EstimatedNvtx) * hgcalSF);
0476         }
0477         if (abs(l1CaloTower.towerEta()) <= 2.1 && abs(l1CaloTower.towerEta()) > 1.8) {
0478           l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() -
0479                                      all_nvtx_to_PU_sub_funcs["hgcalEM"]["er1p8to2p1"].Eval(*EstimatedNvtx) * hgcalSF);
0480         }
0481         if (abs(l1CaloTower.towerEta()) <= 2.4 && abs(l1CaloTower.towerEta()) > 2.1) {
0482           l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() -
0483                                      all_nvtx_to_PU_sub_funcs["hgcalEM"]["er2p1to2p4"].Eval(*EstimatedNvtx) * hgcalSF);
0484         }
0485         if (abs(l1CaloTower.towerEta()) <= 2.7 && abs(l1CaloTower.towerEta()) > 2.4) {
0486           l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() -
0487                                      all_nvtx_to_PU_sub_funcs["hgcalEM"]["er2p4to2p7"].Eval(*EstimatedNvtx) * hgcalSF);
0488         }
0489         if (abs(l1CaloTower.towerEta()) <= 3.1 && abs(l1CaloTower.towerEta()) > 2.7) {
0490           l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() -
0491                                      all_nvtx_to_PU_sub_funcs["hgcalEM"]["er2p7to3p1"].Eval(*EstimatedNvtx) * hgcalSF);
0492         }
0493       }
0494 
0495       // Barrel HCAL eta slices
0496       if (l1CaloTower.hcalTowerEt() > 0. && l1CaloTower.towerIEta() != -98 &&
0497           abs(l1CaloTower.towerEta()) < 2.0)  // abs(eta) < 2 just keeps us out of HF
0498       {
0499         if (abs(l1CaloTower.towerIEta()) <= 3) {
0500           l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() -
0501                                      all_nvtx_to_PU_sub_funcs["hcal"]["er1to3"].Eval(*EstimatedNvtx) * barrelSF);
0502         }
0503         if (abs(l1CaloTower.towerIEta()) <= 6 && abs(l1CaloTower.towerIEta()) >= 4) {
0504           l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() -
0505                                      all_nvtx_to_PU_sub_funcs["hcal"]["er4to6"].Eval(*EstimatedNvtx) * barrelSF);
0506         }
0507         if (abs(l1CaloTower.towerIEta()) <= 9 && abs(l1CaloTower.towerIEta()) >= 7) {
0508           l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() -
0509                                      all_nvtx_to_PU_sub_funcs["hcal"]["er7to9"].Eval(*EstimatedNvtx) * barrelSF);
0510         }
0511         if (abs(l1CaloTower.towerIEta()) <= 12 && abs(l1CaloTower.towerIEta()) >= 10) {
0512           l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() -
0513                                      all_nvtx_to_PU_sub_funcs["hcal"]["er10to12"].Eval(*EstimatedNvtx) * barrelSF);
0514         }
0515         if (abs(l1CaloTower.towerIEta()) <= 15 && abs(l1CaloTower.towerIEta()) >= 13) {
0516           l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() -
0517                                      all_nvtx_to_PU_sub_funcs["hcal"]["er13to15"].Eval(*EstimatedNvtx) * barrelSF);
0518         }
0519         if (abs(l1CaloTower.towerIEta()) <= 18 && abs(l1CaloTower.towerIEta()) >= 16) {
0520           l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() -
0521                                      all_nvtx_to_PU_sub_funcs["hcal"]["er16to18"].Eval(*EstimatedNvtx) * barrelSF);
0522         }
0523       }
0524 
0525       // HGCal Had eta slices
0526       if (l1CaloTower.hcalTowerEt() > 0. && l1CaloTower.towerIEta() == -98) {
0527         if (abs(l1CaloTower.towerEta()) <= 1.8) {
0528           l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() -
0529                                      all_nvtx_to_PU_sub_funcs["hgcalHad"]["er1p4to1p8"].Eval(*EstimatedNvtx) * hgcalSF);
0530         }
0531         if (abs(l1CaloTower.towerEta()) <= 2.1 && abs(l1CaloTower.towerEta()) > 1.8) {
0532           l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() -
0533                                      all_nvtx_to_PU_sub_funcs["hgcalHad"]["er1p8to2p1"].Eval(*EstimatedNvtx) * hgcalSF);
0534         }
0535         if (abs(l1CaloTower.towerEta()) <= 2.4 && abs(l1CaloTower.towerEta()) > 2.1) {
0536           l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() -
0537                                      all_nvtx_to_PU_sub_funcs["hgcalHad"]["er2p1to2p4"].Eval(*EstimatedNvtx) * hgcalSF);
0538         }
0539         if (abs(l1CaloTower.towerEta()) <= 2.7 && abs(l1CaloTower.towerEta()) > 2.4) {
0540           l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() -
0541                                      all_nvtx_to_PU_sub_funcs["hgcalHad"]["er2p4to2p7"].Eval(*EstimatedNvtx) * hgcalSF);
0542         }
0543         if (abs(l1CaloTower.towerEta()) <= 3.1 && abs(l1CaloTower.towerEta()) > 2.7) {
0544           l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() -
0545                                      all_nvtx_to_PU_sub_funcs["hgcalHad"]["er2p7to3p1"].Eval(*EstimatedNvtx) * hgcalSF);
0546         }
0547       }
0548 
0549       // HF eta slices
0550       if (l1CaloTower.hcalTowerEt() > 0. && l1CaloTower.towerIEta() != -98 &&
0551           abs(l1CaloTower.towerEta()) > 2.0)  // abs(eta) > 2 keeps us out of barrel HF
0552       {
0553         if (abs(l1CaloTower.towerIEta()) <= 33 && abs(l1CaloTower.towerIEta()) >= 29) {
0554           l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() -
0555                                      all_nvtx_to_PU_sub_funcs["hf"]["er29to33"].Eval(*EstimatedNvtx) * hfSF);
0556         }
0557         if (abs(l1CaloTower.towerIEta()) <= 37 && abs(l1CaloTower.towerIEta()) >= 34) {
0558           l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() -
0559                                      all_nvtx_to_PU_sub_funcs["hf"]["er34to37"].Eval(*EstimatedNvtx) * hfSF);
0560         }
0561         if (abs(l1CaloTower.towerIEta()) <= 41 && abs(l1CaloTower.towerIEta()) >= 38) {
0562           l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() -
0563                                      all_nvtx_to_PU_sub_funcs["hf"]["er38to41"].Eval(*EstimatedNvtx) * hfSF);
0564         }
0565       }
0566 
0567       // Make sure none are negative
0568       if (l1CaloTower.ecalTowerEt() < 0.)
0569         l1CaloTower.setEcalTowerEt(0.);
0570       if (l1CaloTower.hcalTowerEt() < 0.)
0571         l1CaloTower.setHcalTowerEt(0.);
0572     }
0573   }
0574 
0575   iEvent.put(std::move(EstimatedNvtx), "EstimatedNvtx");
0576   iEvent.put(std::move(L1CaloTowerCalibratedCollection), "L1CaloTowerCalibratedCollection");
0577 }
0578 
0579 DEFINE_FWK_MODULE(L1TowerCalibrator);
[ 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!