L1EGCrystalClusterEmulatorProducer Class Reference

Inheritance diagram for L1EGCrystalClusterEmulatorProducer:

Classes
struct	mycluster
class	SimpleCaloHit

Public Member Functions
	L1EGCrystalClusterEmulatorProducer (const edm::ParameterSet &)
	~L1EGCrystalClusterEmulatorProducer () override
Public Member Functions inherited from edm::stream::EDProducer<>
	EDProducer ()=default
bool	hasAbilityToProduceInBeginLumis () const final
bool	hasAbilityToProduceInBeginProcessBlocks () const final
bool	hasAbilityToProduceInBeginRuns () const final
bool	hasAbilityToProduceInEndLumis () const final
bool	hasAbilityToProduceInEndProcessBlocks () const final
bool	hasAbilityToProduceInEndRuns () const final

Private Member Functions
float	get_calibrate (float uncorr)
bool	passes_iso (float pt, float iso)
bool	passes_looseTkiso (float pt, float iso)
bool	passes_looseTkss (float pt, float ss)
bool	passes_photon (float pt, float pss)
bool	passes_ss (float pt, float ss)
void	produce (edm::Event &, const edm::EventSetup &) override

Private Attributes
l1tp2::ParametricCalibration	calib_
edm::ESHandle< CaloGeometry >	caloGeometry_
edm::ESHandle< CaloTPGTranscoder >	decoder_
const CaloSubdetectorGeometry *	ebGeometry
edm::EDGetTokenT< EcalEBTrigPrimDigiCollection >	ecalTPEBToken_
const CaloSubdetectorGeometry *	hbGeometry
edm::ESHandle< HcalTopology >	hbTopology
edm::EDGetTokenT< edm::SortedCollection< HcalTriggerPrimitiveDigi > >	hcalTPToken_
const HcalTopology *	hcTopology_

Additional Inherited Members
Public Types inherited from edm::stream::EDProducer<>
typedef CacheContexts< T... >	CacheTypes
typedef CacheTypes::GlobalCache	GlobalCache
typedef AbilityChecker< T... >	HasAbility
typedef CacheTypes::LuminosityBlockCache	LuminosityBlockCache
typedef LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >	LuminosityBlockContext
typedef CacheTypes::LuminosityBlockSummaryCache	LuminosityBlockSummaryCache
typedef CacheTypes::RunCache	RunCache
typedef RunContextT< RunCache, GlobalCache >	RunContext
typedef CacheTypes::RunSummaryCache	RunSummaryCache

Detailed Description

Definition at line 226 of file L1EGammaCrystalsEmulatorProducer.cc.

Constructor & Destructor Documentation

L1EGCrystalClusterEmulatorProducer()

L1EGCrystalClusterEmulatorProducer::L1EGCrystalClusterEmulatorProducer ( const edm::ParameterSet &  iConfig )

explicit

Definition at line 339 of file L1EGammaCrystalsEmulatorProducer.cc.

    : ecalTPEBToken_(consumes<EcalEBTrigPrimDigiCollection>(iConfig.getParameter<edm::InputTag>("ecalTPEB"))),
      hcalTPToken_(
          consumes<edm::SortedCollection<HcalTriggerPrimitiveDigi> >(iConfig.getParameter<edm::InputTag>("hcalTP"))),
      calib_(iConfig.getParameter<edm::ParameterSet>("calib")) {
  produces<l1tp2::CaloCrystalClusterCollection>();
  produces<BXVector<l1t::EGamma> >();
  produces<l1tp2::CaloTowerCollection>();
}

~L1EGCrystalClusterEmulatorProducer()

L1EGCrystalClusterEmulatorProducer::~L1EGCrystalClusterEmulatorProducer ( )

override

Definition at line 349 of file L1EGammaCrystalsEmulatorProducer.cc.

{}

Member Function Documentation

get_calibrate()

float L1EGCrystalClusterEmulatorProducer::get_calibrate ( float  uncorr )

private

passes_iso()

bool L1EGCrystalClusterEmulatorProducer::passes_iso ( float  pt, float  iso  )

private

Definition at line 1178 of file L1EGammaCrystalsEmulatorProducer.cc.

                                                                       {
  if (pt < slideIsoPtThreshold) {
    if (!((a0_80 - a1_80 * pt) > iso))
      return false;
  } else {
    if (iso > a0)
      return false;
  }
  return true;
}

References a0, a0_80, a1_80, DiDispStaMuonMonitor_cfi::pt, and slideIsoPtThreshold.

Referenced by produce().

passes_looseTkiso()

bool L1EGCrystalClusterEmulatorProducer::passes_looseTkiso ( float  pt, float  iso  )

private

Definition at line 1189 of file L1EGammaCrystalsEmulatorProducer.cc.

                                                                              {
  return (b0 + b1 * std::exp(-b2 * pt) > iso);
}

References b0, b1, b2, JetChargeProducer_cfi::exp, and DiDispStaMuonMonitor_cfi::pt.

Referenced by produce().

passes_looseTkss()

bool L1EGCrystalClusterEmulatorProducer::passes_looseTkss ( float  pt, float  ss  )

private

Definition at line 1199 of file L1EGammaCrystalsEmulatorProducer.cc.

                                                                            {
  return ((e0_looseTkss - e1_looseTkss * std::exp(-e2_looseTkss * pt)) <= ss);
}

References e0_looseTkss, e1_looseTkss, e2_looseTkss, JetChargeProducer_cfi::exp, DiDispStaMuonMonitor_cfi::pt, and contentValuesCheck::ss.

Referenced by produce().

passes_photon()

bool L1EGCrystalClusterEmulatorProducer::passes_photon ( float  pt, float  pss  )

private

Definition at line 1197 of file L1EGammaCrystalsEmulatorProducer.cc.

{ return (pss > d0 - d1 * pt); }

References d0, d1, and DiDispStaMuonMonitor_cfi::pt.

Referenced by produce().

passes_ss()

bool L1EGCrystalClusterEmulatorProducer::passes_ss ( float  pt, float  ss  )

private

Definition at line 1193 of file L1EGammaCrystalsEmulatorProducer.cc.

                                                                     {
  return ((c0_ss + c1_ss * std::exp(-c2_ss * pt)) <= ss);
}

References c0_ss, c1_ss, c2_ss, JetChargeProducer_cfi::exp, DiDispStaMuonMonitor_cfi::pt, and contentValuesCheck::ss.

Referenced by produce().

produce()

void L1EGCrystalClusterEmulatorProducer::produce ( edm::Event &  iEvent, const edm::EventSetup &  iSetup  )

overrideprivate

Definition at line 351 of file L1EGammaCrystalsEmulatorProducer.cc.

                                                                                              {
  using namespace edm;
 
  edm::Handle<EcalEBTrigPrimDigiCollection> pcalohits;
  iEvent.getByToken(ecalTPEBToken_, pcalohits);
 
  iSetup.get<CaloGeometryRecord>().get(caloGeometry_);
  ebGeometry = caloGeometry_->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
  hbGeometry = caloGeometry_->getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
  iSetup.get<HcalRecNumberingRecord>().get(hbTopology);
  hcTopology_ = hbTopology.product();
  HcalTrigTowerGeometry theTrigTowerGeometry(hcTopology_);
  iSetup.get<CaloTPGRecord>().get(decoder_);
 
  //****************************************************************
  //******************* Get all the hits ***************************
  //****************************************************************
 
  // Get all the ECAL hits
  iEvent.getByToken(ecalTPEBToken_, pcalohits);
  std::vector<SimpleCaloHit> ecalhits;
 
  for (const auto& hit : *pcalohits.product()) {
    if (hit.encodedEt() > 0)  // hit.encodedEt() returns an int corresponding to 2x the crystal Et
    {
      // Et is 10 bit, by keeping the ADC saturation Et at 120 GeV it means that you have to divide by 8
      float et = hit.encodedEt() / 8.;
      if (et < cut_500_MeV)
        continue;  // keep the 500 MeV ET Cut
 
      auto cell = ebGeometry->getGeometry(hit.id());
 
      SimpleCaloHit ehit;
      ehit.setId(hit.id());
      ehit.setPosition(GlobalVector(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z()));
      ehit.setEnergy(et);
      ehit.setPt();
      ecalhits.push_back(ehit);
    }
  }
 
  // Get all the HCAL hits
  std::vector<SimpleCaloHit> hcalhits;
  edm::Handle<edm::SortedCollection<HcalTriggerPrimitiveDigi> > hbhecoll;
  iEvent.getByToken(hcalTPToken_, hbhecoll);
  for (const auto& hit : *hbhecoll.product()) {
    float et = decoder_->hcaletValue(hit.id(), hit.t0());
    if (et <= 0)
      continue;
    if (!(hcTopology_->validHT(hit.id()))) {
      LogError("L1EGCrystalClusterEmulatorProducer")
          << " -- Hcal hit DetID not present in HCAL Geom: " << hit.id() << std::endl;
      throw cms::Exception("L1EGCrystalClusterEmulatorProducer");
      continue;
    }
    const std::vector<HcalDetId>& hcId = theTrigTowerGeometry.detIds(hit.id());
    if (hcId.empty()) {
      LogError("L1EGCrystalClusterEmulatorProducer")
          << "Cannot find any HCalDetId corresponding to " << hit.id() << std::endl;
      throw cms::Exception("L1EGCrystalClusterEmulatorProducer");
      continue;
    }
    if (hcId[0].subdetId() > 1)
      continue;
    GlobalVector hcal_tp_position = GlobalVector(0., 0., 0.);
    for (const auto& hcId_i : hcId) {
      if (hcId_i.subdetId() > 1)
        continue;
      auto cell = hbGeometry->getGeometry(hcId_i);
      if (cell == nullptr)
        continue;
      GlobalVector tmpVector = GlobalVector(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z());
      hcal_tp_position = tmpVector;
      break;
    }
    SimpleCaloHit hhit;
    hhit.setId(hit.id());
    hhit.setIdHcal(hit.id());
    hhit.setPosition(hcal_tp_position);
    hhit.setEnergy(et);
    hhit.setPt();
    hcalhits.push_back(hhit);
  }
 
  //*******************************************************************
  //********************** Do layer 1 *********************************
  //*******************************************************************
 
  // Definition of L1 outputs
  // 36 L1 cards send each 4 links with 17 towers
  float ECAL_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
  float HCAL_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
  int iEta_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
  int iPhi_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
  // 36 L1 cards send each 4 links with 3 clusters
  float energy_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
  // 36 L1 cards send each 4 links with 3 clusters
  int brem_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
  int towerID_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
  int crystalID_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
  int showerShape_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
  int showerShapeLooseTk_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
  int photonShowerShape_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
 
  for (int ii = 0; ii < n_links_card; ++ii) {
    for (int jj = 0; jj < n_towers_per_link; ++jj) {
      for (int ll = 0; ll < n_towers_halfPhi; ++ll) {
        ECAL_tower_L1Card[ii][jj][ll] = 0;
        HCAL_tower_L1Card[ii][jj][ll] = 0;
        iPhi_tower_L1Card[ii][jj][ll] = -999;
        iEta_tower_L1Card[ii][jj][ll] = -999;
      }
    }
  }
  for (int ii = 0; ii < n_links_card; ++ii) {
    for (int jj = 0; jj < n_clusters_link; ++jj) {
      for (int ll = 0; ll < n_towers_halfPhi; ++ll) {
        energy_cluster_L1Card[ii][jj][ll] = 0;
        brem_cluster_L1Card[ii][jj][ll] = 0;
        towerID_cluster_L1Card[ii][jj][ll] = 0;
        crystalID_cluster_L1Card[ii][jj][ll] = 0;
      }
    }
  }
 
  // There is one list of clusters per card. We take the 12 highest pt per card
  vector<mycluster> cluster_list[n_towers_halfPhi];
  // After merging the clusters in different regions of a single L1 card
  vector<mycluster> cluster_list_merged[n_towers_halfPhi];
 
  for (int cc = 0; cc < n_towers_halfPhi; ++cc) {  // Loop over 36 L1 cards
    // Loop over 3x4 etaxphi regions to search for max 5 clusters
    for (int nregion = 0; nregion <= n_clusters_max; ++nregion) {
      int nclusters = 0;
      bool build_cluster = true;
 
      // Continue until 5 clusters have been built or there is no cluster left
      while (nclusters < n_clusters_max && build_cluster) {
        build_cluster = false;
        SimpleCaloHit centerhit;
 
        for (const auto& hit : ecalhits) {
          if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
              getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
              getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
              getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) &&
              // Check that the hit is in the good card
              getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) &&
              getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion &&
              !hit.used() && hit.pt() >= 1.0 && hit.pt() > centerhit.pt())  // 3 towers x 5 crystals
          {
            // Highest hit in good region with pt>1 and not used in any other cluster
            centerhit = hit;
            build_cluster = true;
          }
        }
        if (build_cluster)
          nclusters++;
 
        // Use only the 5 most energetic clusters
        if (build_cluster && nclusters > 0 && nclusters <= n_clusters_max) {
          mycluster mc1;
          mc1.cpt = 0.0;
          mc1.cWeightedEta_ = 0.0;
          mc1.cWeightedPhi_ = 0.0;
          float leftlobe = 0;
          float rightlobe = 0;
          float e5x5 = 0;
          float n5x5 = 0;
          float e2x5_1 = 0;
          float n2x5_1 = 0;
          float e2x5_2 = 0;
          float n2x5_2 = 0;
          float e2x2_1 = 0;
          float n2x2_1 = 0;
          float e2x2_2 = 0;
          float n2x2_2 = 0;
          float e2x2_3 = 0;
          float n2x2_3 = 0;
          float e2x2_4 = 0;
          float n2x2_4 = 0;
          for (auto& hit : ecalhits) {
            if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
                getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
                getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
                getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && hit.pt() > 0 &&
                getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) &&
                getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion) {
              if (abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) > 2 && hit.diphi(centerhit) <= 7) {
                rightlobe += hit.pt();
              }
              if (abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) < -2 && hit.diphi(centerhit) >= -7) {
                leftlobe += hit.pt();
              }
              if (abs(hit.dieta(centerhit)) <= 2 && abs(hit.diphi(centerhit)) <= 2) {
                e5x5 += hit.energy();
                n5x5++;
              }
              if ((hit.dieta(centerhit) == 1 or hit.dieta(centerhit) == 0) &&
                  (hit.diphi(centerhit) == 1 or hit.diphi(centerhit) == 0)) {
                e2x2_1 += hit.energy();
                n2x2_1++;
              }
              if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) &&
                  (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == 1)) {
                e2x2_2 += hit.energy();
                n2x2_2++;
              }
              if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == 1) &&
                  (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == -1)) {
                e2x2_3 += hit.energy();
                n2x2_3++;
              }
              if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) &&
                  (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == -1)) {
                e2x2_4 += hit.energy();
                n2x2_4++;
              }
              if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == 1) && abs(hit.diphi(centerhit)) <= 2) {
                e2x5_1 += hit.energy();
                n2x5_1++;
              }
              if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) && abs(hit.diphi(centerhit)) <= 2) {
                e2x5_2 += hit.energy();
                n2x5_2++;
              }
            }
            if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
                getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
                getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
                getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && !hit.used() && hit.pt() > 0 &&
                abs(hit.dieta(centerhit)) <= 1 && abs(hit.diphi(centerhit)) <= 2 &&
                getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) &&
                getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion) {
              // clusters 3x5 in etaxphi using only the hits in the corresponding card and in the corresponding 3x4 region
              hit.setUsed(true);
              mc1.cpt += hit.pt();
              mc1.cWeightedEta_ += float(hit.pt()) * float(hit.position().eta());
              mc1.cWeightedPhi_ = mc1.cWeightedPhi_ + (float(hit.pt()) * float(hit.position().phi()));
            }
          }
          if (do_brem && (rightlobe > 0.10 * mc1.cpt or leftlobe > 0.10 * mc1.cpt)) {
            for (auto& hit : ecalhits) {
              if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
                  getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
                  getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
                  getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && hit.pt() > 0 &&
                  getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) &&
                  getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion &&
                  !hit.used()) {
                if (rightlobe > 0.10 * mc1.cpt && (leftlobe < 0.10 * mc1.cpt or rightlobe > leftlobe) &&
                    abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) > 2 && hit.diphi(centerhit) <= 7) {
                  mc1.cpt += hit.pt();
                  hit.setUsed(true);
                  mc1.cbrem_ = 1;
                }
                if (leftlobe > 0.10 * mc1.cpt && (rightlobe < 0.10 * mc1.cpt or leftlobe >= rightlobe) &&
                    abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) < -2 && hit.diphi(centerhit) >= -7) {
                  mc1.cpt += hit.pt();
                  hit.setUsed(true);
                  mc1.cbrem_ = 1;
                }
              }
            }
          }
          mc1.c5x5_ = e5x5;
          mc1.c2x5_ = max(e2x5_1, e2x5_2);
          mc1.c2x2_ = e2x2_1;
          if (e2x2_2 > mc1.c2x2_)
            mc1.c2x2_ = e2x2_2;
          if (e2x2_3 > mc1.c2x2_)
            mc1.c2x2_ = e2x2_3;
          if (e2x2_4 > mc1.c2x2_)
            mc1.c2x2_ = e2x2_4;
          mc1.cWeightedEta_ = mc1.cWeightedEta_ / mc1.cpt;
          mc1.cWeightedPhi_ = mc1.cWeightedPhi_ / mc1.cpt;
          mc1.ceta_ = getCrystal_etaID(centerhit.position().eta());
          mc1.cphi_ = getCrystal_phiID(centerhit.position().phi());
          mc1.crawphi_ = centerhit.position().phi();
          mc1.craweta_ = centerhit.position().eta();
          cluster_list[cc].push_back(mc1);
        }  // End if 5 clusters per region
      }    // End while to find the 5 clusters
    }      // End loop over regions to search for clusters
    std::sort(begin(cluster_list[cc]), end(cluster_list[cc]), [](mycluster a, mycluster b) { return a.cpt > b.cpt; });
 
    // Merge clusters from different regions
    for (unsigned int jj = 0; jj < unsigned(cluster_list[cc].size()); ++jj) {
      for (unsigned int kk = jj + 1; kk < unsigned(cluster_list[cc].size()); ++kk) {
        if (std::abs(cluster_list[cc][jj].ceta_ - cluster_list[cc][kk].ceta_) < 2 &&
            std::abs(cluster_list[cc][jj].cphi_ - cluster_list[cc][kk].cphi_) < 2) {  //Diagonale + exact neighbors
          if (cluster_list[cc][kk].cpt > cluster_list[cc][jj].cpt) {
            cluster_list[cc][kk].cpt += cluster_list[cc][jj].cpt;
            cluster_list[cc][kk].c5x5_ += cluster_list[cc][jj].c5x5_;
            cluster_list[cc][kk].c2x5_ += cluster_list[cc][jj].c2x5_;
            cluster_list[cc][jj].cpt = 0;
            cluster_list[cc][jj].c5x5_ = 0;
            cluster_list[cc][jj].c2x5_ = 0;
            cluster_list[cc][jj].c2x2_ = 0;
          } else {
            cluster_list[cc][jj].cpt += cluster_list[cc][kk].cpt;
            cluster_list[cc][jj].c5x5_ += cluster_list[cc][kk].c5x5_;
            cluster_list[cc][jj].c2x5_ += cluster_list[cc][kk].c2x5_;
            cluster_list[cc][kk].cpt = 0;
            cluster_list[cc][kk].c2x2_ = 0;
            cluster_list[cc][kk].c2x5_ = 0;
            cluster_list[cc][kk].c5x5_ = 0;
          }
        }
      }
      if (cluster_list[cc][jj].cpt > 0) {
        cluster_list[cc][jj].cpt =
            cluster_list[cc][jj].cpt *
            calib_(cluster_list[cc][jj].cpt,
                   std::abs(cluster_list[cc][jj].craweta_));  //Mark's calibration as a function of eta and pt
        cluster_list_merged[cc].push_back(cluster_list[cc][jj]);
      }
    }
    std::sort(begin(cluster_list_merged[cc]), end(cluster_list_merged[cc]), [](mycluster a, mycluster b) {
      return a.cpt > b.cpt;
    });
 
    // Fill cluster information in the arrays. We keep max 12 clusters (distributed between 4 links)
    for (unsigned int jj = 0; jj < unsigned(cluster_list_merged[cc].size()) && jj < n_clusters_4link; ++jj) {
      crystalID_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] =
          getCrystalIDInTower(cluster_list_merged[cc][jj].ceta_, cluster_list_merged[cc][jj].cphi_);
      towerID_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] =
          getTowerID(cluster_list_merged[cc][jj].ceta_, cluster_list_merged[cc][jj].cphi_);
      energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = cluster_list_merged[cc][jj].cpt;
      brem_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = cluster_list_merged[cc][jj].cbrem_;
      if (passes_ss(cluster_list_merged[cc][jj].cpt,
                    cluster_list_merged[cc][jj].c2x5_ / cluster_list_merged[cc][jj].c5x5_))
        showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1;
      else
        showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0;
      if (passes_looseTkss(cluster_list_merged[cc][jj].cpt,
                           cluster_list_merged[cc][jj].c2x5_ / cluster_list_merged[cc][jj].c5x5_))
        showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1;
      else
        showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0;
      if (passes_photon(cluster_list_merged[cc][jj].cpt,
                        cluster_list_merged[cc][jj].c2x2_ / cluster_list_merged[cc][jj].c2x5_))
        photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1;
      else
        photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0;
    }
 
    // Loop over calo ecal hits to get the ECAL towers. Take only hits that have not been used to make clusters
    for (const auto& hit : ecalhits) {
      if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
          getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
          getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
          getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) &&
          !hit.used()) {                             // Take all the hits inside the card that have not been used yet
        for (int jj = 0; jj < n_links_card; ++jj) {  // loop over 4 links per card
          if ((getCrystal_phiID(hit.position().phi()) / n_crystals_towerPhi) % 4 == jj) {  // Go to ID tower modulo 4
            for (int ii = 0; ii < n_towers_per_link; ++ii) {
              //Apply Mark's calibration at the same time (row of the lowest pT, as a function of eta)
              if ((getCrystal_etaID(hit.position().eta()) / n_crystals_towerEta) % n_towers_per_link == ii) {
                ECAL_tower_L1Card[jj][ii][cc] += hit.pt() * calib_(0, std::abs(hit.position().eta()));
                iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(hit.position().eta());  //hit.id().ieta();
                iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(hit.position().phi());  //hit.id().iphi();
              }
            }  // end of loop over eta towers
          }
        }  // end of loop over phi links
        // Make sure towers with 0 ET are initialized with proper iEta, iPhi coordinates
        static constexpr float tower_width = 0.0873;
        for (int jj = 0; jj < n_links_card; ++jj) {
          for (int ii = 0; ii < n_towers_per_link; ++ii) {
            float phi = getPhiMin_card(cc) * tower_width / n_crystals_towerPhi - M_PI + (jj + 0.5) * tower_width;
            float eta = getEtaMin_card(cc) * tower_width / n_crystals_towerEta - n_towers_per_link * tower_width +
                        (ii + 0.5) * tower_width;
            iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(eta);
            iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(phi);
          }
        }
 
      }  // end of check if inside card
    }    // end of loop over hits to build towers
 
    // Loop over hcal hits to get the HCAL towers.
    for (const auto& hit : hcalhits) {
      if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
          getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
          getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
          getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && hit.pt() > 0) {
        for (int jj = 0; jj < n_links_card; ++jj) {
          if ((getCrystal_phiID(hit.position().phi()) / n_crystals_towerPhi) % n_links_card == jj) {
            for (int ii = 0; ii < n_towers_per_link; ++ii) {
              if ((getCrystal_etaID(hit.position().eta()) / n_crystals_towerEta) % n_towers_per_link == ii) {
                HCAL_tower_L1Card[jj][ii][cc] += hit.pt();
                iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(hit.position().eta());  //hit.id().ieta();
                iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(hit.position().phi());  //hit.id().iphi();
              }
            }  // end of loop over eta towers
          }
        }  // end of loop over phi links
      }    // end of check if inside card
    }      // end of loop over hits to build towers
 
    // Give back energy of not used clusters to the towers (if there are more than 12 clusters)
    for (unsigned int kk = n_clusters_4link; kk < cluster_list_merged[cc].size(); ++kk) {
      if (cluster_list_merged[cc][kk].cpt > 0) {
        ECAL_tower_L1Card[getTower_phiID(cluster_list_merged[cc][kk].cphi_)]
                         [getTower_etaID(cluster_list_merged[cc][kk].ceta_)][cc] += cluster_list_merged[cc][kk].cpt;
      }
    }
  }  //End of loop over cards
 
  //*********************************************************
  //******************** Do layer 2 *************************
  //*********************************************************
 
  // Definition of L2 outputs
  float HCAL_tower_L2Card[n_links_GCTcard][n_towers_per_link]
                         [n_GCTcards];  // 3 L2 cards send each 48 links with 17 towers
  float ECAL_tower_L2Card[n_links_GCTcard][n_towers_per_link][n_GCTcards];
  int iEta_tower_L2Card[n_links_GCTcard][n_towers_per_link][n_GCTcards];
  int iPhi_tower_L2Card[n_links_GCTcard][n_towers_per_link][n_GCTcards];
  float energy_cluster_L2Card[n_links_GCTcard][n_clusters_per_link]
                             [n_GCTcards];  // 3 L2 cards send each 48 links with 2 clusters
  float brem_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
  int towerID_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
  int crystalID_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
  float isolation_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
  float HE_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
  int showerShape_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
  int showerShapeLooseTk_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
  int photonShowerShape_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
 
  for (int ii = 0; ii < n_links_GCTcard; ++ii) {
    for (int jj = 0; jj < n_towers_per_link; ++jj) {
      for (int ll = 0; ll < n_GCTcards; ++ll) {
        HCAL_tower_L2Card[ii][jj][ll] = 0;
        ECAL_tower_L2Card[ii][jj][ll] = 0;
        iEta_tower_L2Card[ii][jj][ll] = 0;
        iPhi_tower_L2Card[ii][jj][ll] = 0;
      }
    }
  }
  for (int ii = 0; ii < n_links_GCTcard; ++ii) {
    for (int jj = 0; jj < n_clusters_per_link; ++jj) {
      for (int ll = 0; ll < n_GCTcards; ++ll) {
        energy_cluster_L2Card[ii][jj][ll] = 0;
        brem_cluster_L2Card[ii][jj][ll] = 0;
        towerID_cluster_L2Card[ii][jj][ll] = 0;
        crystalID_cluster_L2Card[ii][jj][ll] = 0;
        isolation_cluster_L2Card[ii][jj][ll] = 0;
        HE_cluster_L2Card[ii][jj][ll] = 0;
        photonShowerShape_cluster_L2Card[ii][jj][ll] = 0;
        showerShape_cluster_L2Card[ii][jj][ll] = 0;
        showerShapeLooseTk_cluster_L2Card[ii][jj][ll] = 0;
      }
    }
  }
 
  // There is one list of clusters per equivalent of L1 card. We take the 8 highest pt.
  vector<mycluster> cluster_list_L2[n_towers_halfPhi];
 
  // Merge clusters on the phi edges
  for (int ii = 0; ii < n_borders_phi; ++ii) {  // 18 borders in phi
    for (int jj = 0; jj < n_eta_bins; ++jj) {   // 2 eta bins
      int card_left = 2 * ii + jj;
      int card_right = 2 * ii + jj + 2;
      if (card_right > 35)
        card_right = card_right - n_towers_halfPhi;
      for (int kk = 0; kk < n_clusters_4link; ++kk) {  // 12 clusters in the first card. We check the right side
        if (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 50 &&
            crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 19 &&
            energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 0) {
          for (int ll = 0; ll < n_clusters_4link; ++ll) {  // We check the 12 clusters in the card on the right
            if (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < n_towers_per_link &&
                crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < 5 &&
                std::abs(
                    5 * (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right]) % n_towers_per_link +
                    crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] % 5 -
                    5 * (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left]) % n_towers_per_link -
                    crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] % 5) < 2) {
              if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] >
                  energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right]) {
                energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] +=
                    energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right];
                energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] = 0;
              }  // The least energetic cluster is merged to the most energetic one
              else {
                energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] +=
                    energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left];
                energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] = 0;
              }
            }
          }
        }
      }
    }
  }
 
  // Bremsstrahlung corrections. Merge clusters on the phi edges depending on pT (pt more than 10 percent, dphi leq 5, deta leq 1)
  if (do_brem) {
    for (int ii = 0; ii < n_borders_phi; ++ii) {  // 18 borders in phi
      for (int jj = 0; jj < n_eta_bins; ++jj) {   // 2 eta bins
        int card_left = 2 * ii + jj;
        int card_right = 2 * ii + jj + 2;
        if (card_right > 35)
          card_right = card_right - n_towers_halfPhi;
        // 12 clusters in the first card. We check the right side crystalID_cluster_L1Card[kk%4][kk/4][card_left]
        for (int kk = 0; kk < n_clusters_4link; ++kk) {
          // if the tower is at the edge there might be brem corrections, whatever the crystal ID
          if (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 50 &&
              energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 0) {
            for (int ll = 0; ll < n_clusters_4link; ++ll) {  // We check the 12 clusters in the card on the right
              //Distance of 1 max in eta
              if (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < n_towers_per_link &&
                  std::abs(5 * (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right]) %
                               n_towers_per_link +
                           crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] % 5 -
                           5 * (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left]) %
                               n_towers_per_link -
                           crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] % 5) <= 1) {
                //Distance of 5 max in phi
                if (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < n_towers_per_link &&
                    std::abs(5 + crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] / 5 -
                             (crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] / 5)) <= 5) {
                  if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] >
                          energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] &&
                      energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] >
                          0.10 * energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left]) {
                    energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] +=
                        energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right];
                    energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] = 0;
                    brem_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] = 1;
                  }
                  // The least energetic cluster is merged to the most energetic one, if its pt is at least ten percent
                  else if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_right] >=
                               energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_left] &&
                           energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_left] >
                               0.10 * energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_right]) {
                    energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] +=
                        energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left];
                    energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] = 0;
                    brem_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] = 1;
                  }
                }  //max distance eta
              }    //max distance phi
            }      //max distance phi
          }
        }
      }
    }
  }
 
  // Merge clusters on the eta edges
  for (int ii = 0; ii < n_borders_eta; ++ii) {  // 18 borders in eta
    int card_bottom = 2 * ii;
    int card_top = 2 * ii + 1;
    for (int kk = 0; kk < n_clusters_4link; ++kk) {  // 12 clusters in the first card. We check the top side
      if (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] % n_towers_per_link == 16 &&
          crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] % 5 == n_links_card &&
          energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] >
              0) {                                       // If there is one cluster on the right side of the first card
        for (int ll = 0; ll < n_clusters_4link; ++ll) {  // We check the card on the right
          if (std::abs(
                  5 * (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] / n_towers_per_link) +
                  crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] / 5 -
                  5 * (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] / n_towers_per_link) -
                  crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] / 5) < 2) {
            if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] >
                energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_bottom]) {
              energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] +=
                  energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top];
              energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] = 0;
            } else {
              energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] +=
                  energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom];
              energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] = 0;
            }
          }
        }
      }
    }
  }
 
  // Regroup the new clusters per equivalent of L1 card geometry
  for (int ii = 0; ii < n_towers_halfPhi; ++ii) {
    for (int jj = 0; jj < n_clusters_4link; ++jj) {
      if (energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii] > 0) {
        mycluster mc1;
        mc1.cpt = energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
        mc1.cbrem_ = brem_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
        mc1.ctowerid_ = towerID_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
        mc1.ccrystalid_ = crystalID_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
        mc1.cshowershape_ = showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
        mc1.cshowershapeloosetk_ = showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
        mc1.cphotonshowershape_ = photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
        cluster_list_L2[ii].push_back(mc1);
      }
    }
    std::sort(
        begin(cluster_list_L2[ii]), end(cluster_list_L2[ii]), [](mycluster a, mycluster b) { return a.cpt > b.cpt; });
  }
 
  // If there are more than 8 clusters per equivalent of L1 card we need to put them back in the towers
  for (int ii = 0; ii < n_towers_halfPhi; ++ii) {
    for (unsigned int jj = 8; jj < n_clusters_4link && jj < cluster_list_L2[ii].size(); ++jj) {
      if (cluster_list_L2[ii][jj].cpt > 0) {
        ECAL_tower_L1Card[cluster_list_L2[ii][jj].ctowerid_ / n_towers_per_link]
                         [cluster_list_L2[ii][jj].ctowerid_ % n_towers_per_link][ii] += cluster_list_L2[ii][jj].cpt;
        cluster_list_L2[ii][jj].cpt = 0;
        cluster_list_L2[ii][jj].ctowerid_ = 0;
        cluster_list_L2[ii][jj].ccrystalid_ = 0;
      }
    }
  }
 
  // Compute isolation (7*7 ECAL towers) and HCAL energy (5x5 HCAL towers)
  for (int ii = 0; ii < n_towers_halfPhi; ++ii) {  // Loop over the new cluster list (stored in 36x8 format)
    for (unsigned int jj = 0; jj < 8 && jj < cluster_list_L2[ii].size(); ++jj) {
      int cluster_etaOfTower_fullDetector = get_towerEta_fromCardTowerInCard(ii, cluster_list_L2[ii][jj].ctowerid_);
      int cluster_phiOfTower_fullDetector = get_towerPhi_fromCardTowerInCard(ii, cluster_list_L2[ii][jj].ctowerid_);
      float hcal_nrj = 0.0;
      float isolation = 0.0;
      int ntowers = 0;
      for (int iii = 0; iii < n_towers_halfPhi; ++iii) {  // The clusters have to be added to the isolation
        for (unsigned int jjj = 0; jjj < 8 && jjj < cluster_list_L2[iii].size(); ++jjj) {
          if (!(iii == ii && jjj == jj)) {
            int cluster2_eta = get_towerEta_fromCardTowerInCard(iii, cluster_list_L2[iii][jjj].ctowerid_);
            int cluster2_phi = get_towerPhi_fromCardTowerInCard(iii, cluster_list_L2[iii][jjj].ctowerid_);
            if (abs(cluster2_eta - cluster_etaOfTower_fullDetector) <= 2 &&
                (abs(cluster2_phi - cluster_phiOfTower_fullDetector) <= 2 or
                 abs(cluster2_phi - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) {
              isolation += cluster_list_L2[iii][jjj].cpt;
            }
          }
        }
      }
      for (int kk = 0; kk < n_towers_halfPhi; ++kk) {       // 36 cards
        for (int ll = 0; ll < n_links_card; ++ll) {         // 4 links per card
          for (int mm = 0; mm < n_towers_per_link; ++mm) {  // 17 towers per link
            int etaOftower_fullDetector = get_towerEta_fromCardLinkTower(kk, ll, mm);
            int phiOftower_fullDetector = get_towerPhi_fromCardLinkTower(kk, ll, mm);
            // First do ECAL
            // The towers are within 3. Needs to stitch the two phi sides together
            if (abs(etaOftower_fullDetector - cluster_etaOfTower_fullDetector) <= 2 &&
                (abs(phiOftower_fullDetector - cluster_phiOfTower_fullDetector) <= 2 or
                 abs(phiOftower_fullDetector - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) {
              // Remove the column outside of the L2 card:  values (0,71), (23,26), (24,21), (47,50), (48,50), (71,2)
              if (!((cluster_phiOfTower_fullDetector == 0 && phiOftower_fullDetector == 71) or
                    (cluster_phiOfTower_fullDetector == 23 && phiOftower_fullDetector == 26) or
                    (cluster_phiOfTower_fullDetector == 24 && phiOftower_fullDetector == 21) or
                    (cluster_phiOfTower_fullDetector == 47 && phiOftower_fullDetector == 50) or
                    (cluster_phiOfTower_fullDetector == 48 && phiOftower_fullDetector == 45) or
                    (cluster_phiOfTower_fullDetector == 71 && phiOftower_fullDetector == 2))) {
                isolation += ECAL_tower_L1Card[ll][mm][kk];
                ntowers++;
              }
            }
            // Now do HCAL
            // The towers are within 2. Needs to stitch the two phi sides together
            if (abs(etaOftower_fullDetector - cluster_etaOfTower_fullDetector) <= 2 &&
                (abs(phiOftower_fullDetector - cluster_phiOfTower_fullDetector) <= 2 or
                 abs(phiOftower_fullDetector - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) {
              hcal_nrj += HCAL_tower_L1Card[ll][mm][kk];
            }
          }
        }
      }
      cluster_list_L2[ii][jj].ciso_ = ((isolation) * (25.0 / ntowers)) / cluster_list_L2[ii][jj].cpt;
      cluster_list_L2[ii][jj].chovere_ = hcal_nrj / cluster_list_L2[ii][jj].cpt;
    }
  }
 
  //Reformat the information inside the 3 L2 cards
  //First let's fill the towers
  for (int ii = 0; ii < n_links_GCTcard; ++ii) {
    for (int jj = 0; jj < n_towers_per_link; ++jj) {
      for (int ll = 0; ll < 3; ++ll) {
        ECAL_tower_L2Card[ii][jj][ll] =
            ECAL_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
        HCAL_tower_L2Card[ii][jj][ll] =
            HCAL_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
        iEta_tower_L2Card[ii][jj][ll] =
            iEta_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
        iPhi_tower_L2Card[ii][jj][ll] =
            iPhi_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
      }
    }
  }
 
  //Second let's fill the clusters
  for (int ii = 0; ii < n_towers_halfPhi; ++ii) {  // The cluster list is still in the L1 like geometry
    for (unsigned int jj = 0; jj < unsigned(cluster_list_L2[ii].size()) && jj < n_clusters_4link; ++jj) {
      crystalID_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
                              [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ccrystalid_;
      towerID_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
                            [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ctowerid_;
      energy_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
                           [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cpt;
      brem_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
                         [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cbrem_;
      isolation_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
                              [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ciso_;
      HE_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
                       [ii / n_clusters_4link] = cluster_list_L2[ii][jj].chovere_;
      showerShape_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
                                [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cshowershape_;
      showerShapeLooseTk_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
                                       [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cshowershapeloosetk_;
      photonShowerShape_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
                                      [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cphotonshowershape_;
    }
  }
 
  auto L1EGXtalClusters = std::make_unique<l1tp2::CaloCrystalClusterCollection>();
  auto L1EGammas = std::make_unique<l1t::EGammaBxCollection>();
  auto L1CaloTowers = std::make_unique<l1tp2::CaloTowerCollection>();
 
  // Fill the cluster collection and towers as well
  for (int ii = 0; ii < n_links_GCTcard; ++ii) {  // 48 links
    for (int ll = 0; ll < n_GCTcards; ++ll) {     // 3 cards
      // For looping over the Towers a few lines below
      for (int jj = 0; jj < 2; ++jj) {  // 2 L1EGs
        if (energy_cluster_L2Card[ii][jj][ll] > 0.45) {
          reco::Candidate::PolarLorentzVector p4calibrated(
              energy_cluster_L2Card[ii][jj][ll],
              getEta_fromL2LinkCardTowerCrystal(
                  ii, ll, towerID_cluster_L2Card[ii][jj][ll], crystalID_cluster_L2Card[ii][jj][ll]),
              getPhi_fromL2LinkCardTowerCrystal(
                  ii, ll, towerID_cluster_L2Card[ii][jj][ll], crystalID_cluster_L2Card[ii][jj][ll]),
              0.);
          SimpleCaloHit centerhit;
          bool is_iso = passes_iso(energy_cluster_L2Card[ii][jj][ll], isolation_cluster_L2Card[ii][jj][ll]);
          bool is_looseTkiso =
              passes_looseTkiso(energy_cluster_L2Card[ii][jj][ll], isolation_cluster_L2Card[ii][jj][ll]);
          bool is_ss = (showerShape_cluster_L2Card[ii][jj][ll] == 1);
          bool is_photon = (photonShowerShape_cluster_L2Card[ii][jj][ll] == 1) && is_ss && is_iso;
          bool is_looseTkss = (showerShapeLooseTk_cluster_L2Card[ii][jj][ll] == 1);
          // All the ID set to Standalone WP! Some dummy values for non calculated variables
          l1tp2::CaloCrystalCluster cluster(p4calibrated,
                                            energy_cluster_L2Card[ii][jj][ll],
                                            HE_cluster_L2Card[ii][jj][ll],
                                            isolation_cluster_L2Card[ii][jj][ll],
                                            centerhit.id(),
                                            -1000,
                                            float(brem_cluster_L2Card[ii][jj][ll]),
                                            -1000,
                                            -1000,
                                            energy_cluster_L2Card[ii][jj][ll],
                                            -1,
                                            is_iso&& is_ss,
                                            is_iso&& is_ss,
                                            is_photon,
                                            is_iso&& is_ss,
                                            is_looseTkiso&& is_looseTkss,
                                            is_iso&& is_ss);
          // Experimental parameters, don't want to bother with hardcoding them in data format
          std::map<std::string, float> params;
          params["standaloneWP_showerShape"] = is_ss;
          params["standaloneWP_isolation"] = is_iso;
          params["trkMatchWP_showerShape"] = is_looseTkss;
          params["trkMatchWP_isolation"] = is_looseTkiso;
          params["TTiEta"] = getToweriEta_fromAbsoluteID(getTower_absoluteEtaID(p4calibrated.eta()));
          params["TTiPhi"] = getToweriPhi_fromAbsoluteID(getTower_absolutePhiID(p4calibrated.phi()));
          cluster.setExperimentalParams(params);
          L1EGXtalClusters->push_back(cluster);
 
          int standaloneWP = (int)(is_iso && is_ss);
          int looseL1TkMatchWP = (int)(is_looseTkiso && is_looseTkss);
          int photonWP = (int)(is_photon);
          int quality =
              (standaloneWP * std::pow(2, 0)) + (looseL1TkMatchWP * std::pow(2, 1)) + (photonWP * std::pow(2, 2));
          L1EGammas->push_back(
              0, l1t::EGamma(p4calibrated, p4calibrated.pt(), p4calibrated.eta(), p4calibrated.phi(), quality, 1));
        }
      }
    }
  }
 
  for (int ii = 0; ii < n_links_GCTcard; ++ii) {  // 48 links
    for (int ll = 0; ll < n_GCTcards; ++ll) {     // 3 cards
      // Fill the tower collection
      for (int jj = 0; jj < n_towers_per_link; ++jj) {  // 17 TTs
        l1tp2::CaloTower l1CaloTower;
        l1CaloTower.setEcalTowerEt(ECAL_tower_L2Card[ii][jj][ll]);
        l1CaloTower.setHcalTowerEt(HCAL_tower_L2Card[ii][jj][ll]);
        l1CaloTower.setTowerIEta(getToweriEta_fromAbsoluteID(iEta_tower_L2Card[ii][jj][ll]));
        l1CaloTower.setTowerIPhi(getToweriPhi_fromAbsoluteID(iPhi_tower_L2Card[ii][jj][ll]));
        l1CaloTower.setTowerEta(getTowerEta_fromAbsoluteID(iEta_tower_L2Card[ii][jj][ll]));
        l1CaloTower.setTowerPhi(getTowerPhi_fromAbsoluteID(iPhi_tower_L2Card[ii][jj][ll]));
        // Some towers have incorrect eta/phi because that wasn't initialized in certain cases, fix these
        static float constexpr towerEtaUpperUnitialized = -80;
        static float constexpr towerPhiUpperUnitialized = -90;
        if (l1CaloTower.towerEta() < towerEtaUpperUnitialized && l1CaloTower.towerPhi() < towerPhiUpperUnitialized) {
          l1CaloTower.setTowerEta(l1t::CaloTools::towerEta(l1CaloTower.towerIEta()));
          l1CaloTower.setTowerPhi(l1t::CaloTools::towerPhi(l1CaloTower.towerIEta(), l1CaloTower.towerIPhi()));
        }
        l1CaloTower.setIsBarrel(true);
 
        // Add L1EGs if they match in iEta / iPhi
        // L1EGs are already pT ordered, we will take the ID info for the leading one, but pT as the sum
        for (const auto& l1eg : *L1EGXtalClusters) {
          if (l1eg.experimentalParam("TTiEta") != l1CaloTower.towerIEta())
            continue;
          if (l1eg.experimentalParam("TTiPhi") != l1CaloTower.towerIPhi())
            continue;
 
          int n_L1eg = l1CaloTower.nL1eg();
          l1CaloTower.setNL1eg(n_L1eg++);
          l1CaloTower.setL1egTowerEt(l1CaloTower.l1egTowerEt() + l1eg.pt());
          // Don't record L1EG quality info for subleading L1EG
          if (l1CaloTower.nL1eg() > 1)
            continue;
          l1CaloTower.setL1egTrkSS(l1eg.experimentalParam("trkMatchWP_showerShape"));
          l1CaloTower.setL1egTrkIso(l1eg.experimentalParam("trkMatchWP_isolation"));
          l1CaloTower.setL1egStandaloneSS(l1eg.experimentalParam("standaloneWP_showerShape"));
          l1CaloTower.setL1egStandaloneIso(l1eg.experimentalParam("standaloneWP_isolation"));
        }
 
        L1CaloTowers->push_back(l1CaloTower);
      }
    }
  }
 
  iEvent.put(std::move(L1EGXtalClusters));
  iEvent.put(std::move(L1EGammas));
  iEvent.put(std::move(L1CaloTowers));
}

References a, funct::abs(), b, L1EGCrystalClusterEmulatorProducer::mycluster::c2x2_, L1EGCrystalClusterEmulatorProducer::mycluster::c2x5_, L1EGCrystalClusterEmulatorProducer::mycluster::c5x5_, calib_, caloGeometry_, L1EGCrystalClusterEmulatorProducer::mycluster::cbrem_, L1EGCrystalClusterEmulatorProducer::mycluster::ccrystalid_, L1EGCrystalClusterEmulatorProducer::mycluster::ceta_, convert_L2toL1_card(), convert_L2toL1_link(), convert_L2toL1_tower(), L1EGCrystalClusterEmulatorProducer::mycluster::cphi_, L1EGCrystalClusterEmulatorProducer::mycluster::cphotonshowershape_, L1EGCrystalClusterEmulatorProducer::mycluster::cpt, L1EGCrystalClusterEmulatorProducer::mycluster::craweta_, L1EGCrystalClusterEmulatorProducer::mycluster::crawphi_, L1EGCrystalClusterEmulatorProducer::mycluster::cshowershape_, L1EGCrystalClusterEmulatorProducer::mycluster::cshowershapeloosetk_, L1EGCrystalClusterEmulatorProducer::mycluster::ctowerid_, cut_500_MeV, L1EGCrystalClusterEmulatorProducer::mycluster::cWeightedEta_, L1EGCrystalClusterEmulatorProducer::mycluster::cWeightedPhi_, decoder_, HcalTrigTowerGeometry::detIds(), do_brem, ebGeometry, DetId::Ecal, EcalBarrel, ecalTPEBToken_, mps_fire::end, EgHLTOffHistBins_cfi::et, PVValHelper::eta, PV3DBase< T, PVType, FrameType >::eta(), Exception, dqmMemoryStats::float, edm::EventSetup::get(), get, get_towerEta_fromCardLinkTower(), get_towerEta_fromCardTowerInCard(), get_towerPhi_fromCardLinkTower(), get_towerPhi_fromCardTowerInCard(), getCrystal_etaID(), getCrystal_phiID(), getCrystalIDInTower(), getEta_fromL2LinkCardTowerCrystal(), getEtaMax_card(), getEtaMin_card(), CaloSubdetectorGeometry::getGeometry(), getPhi_fromL2LinkCardTowerCrystal(), getPhiMax_card(), getPhiMin_card(), CaloGeometry::getSubdetectorGeometry(), getTower_absoluteEtaID(), getTower_absolutePhiID(), getTower_etaID(), getTower_phiID(), getTowerEta_fromAbsoluteID(), getTowerID(), getToweriEta_fromAbsoluteID(), getToweriPhi_fromAbsoluteID(), getTowerPhi_fromAbsoluteID(), hbGeometry, hbTopology, DetId::Hcal, HcalBarrel, CaloTPGTranscoder::hcaletValue(), hcalTPToken_, hcTopology_, hit::id, L1EGCrystalClusterEmulatorProducer::SimpleCaloHit::id(), iEvent, cuy::ii, createfilelist::int, findQualityFiles::jj, GetRecoTauVFromDQM_MC_cff::kk, l1tp2::CaloTower::l1egTowerEt(), M_PI, SiStripPI::max, eostools::move(), n_borders_eta, n_borders_phi, n_clusters_4link, n_clusters_link, n_clusters_max, n_clusters_per_link, n_crystals_3towers, n_crystals_towerEta, n_crystals_towerPhi, n_eta_bins, n_GCTcards, n_links_card, n_links_GCTcard, n_towers_halfPhi, n_towers_per_link, n_towers_Phi, l1tp2::CaloTower::nL1eg(), or, CalibrationSummaryClient_cfi::params, passes_iso(), passes_looseTkiso(), passes_looseTkss(), passes_photon(), passes_ss(), phi, PV3DBase< T, PVType, FrameType >::phi(), L1EGCrystalClusterEmulatorProducer::SimpleCaloHit::position(), funct::pow(), edm::Handle< T >::product(), edm::ESHandle< T >::product(), L1EGCrystalClusterEmulatorProducer::SimpleCaloHit::pt(), qcdUeDQM_cfi::quality, l1tp2::CaloTower::setEcalTowerEt(), L1EGCrystalClusterEmulatorProducer::SimpleCaloHit::setEnergy(), l1tp2::CaloTower::setHcalTowerEt(), L1EGCrystalClusterEmulatorProducer::SimpleCaloHit::setId(), L1EGCrystalClusterEmulatorProducer::SimpleCaloHit::setIdHcal(), l1tp2::CaloTower::setIsBarrel(), l1tp2::CaloTower::setL1egStandaloneIso(), l1tp2::CaloTower::setL1egStandaloneSS(), l1tp2::CaloTower::setL1egTowerEt(), l1tp2::CaloTower::setL1egTrkIso(), l1tp2::CaloTower::setL1egTrkSS(), l1tp2::CaloTower::setNL1eg(), L1EGCrystalClusterEmulatorProducer::SimpleCaloHit::setPosition(), L1EGCrystalClusterEmulatorProducer::SimpleCaloHit::setPt(), l1tp2::CaloTower::setTowerEta(), l1tp2::CaloTower::setTowerIEta(), l1tp2::CaloTower::setTowerIPhi(), l1tp2::CaloTower::setTowerPhi(), findQualityFiles::size, l1tp2::CaloTower::towerEta(), l1t::CaloTools::towerEta(), l1tp2::CaloTower::towerIEta(), l1tp2::CaloTower::towerIPhi(), l1tp2::CaloTower::towerPhi(), l1t::CaloTools::towerPhi(), and HcalTopology::validHT().

Member Data Documentation

calib_

l1tp2::ParametricCalibration L1EGCrystalClusterEmulatorProducer::calib_

private

Definition at line 244 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

caloGeometry_

edm::ESHandle<CaloGeometry> L1EGCrystalClusterEmulatorProducer::caloGeometry_

private

Definition at line 246 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

decoder_

edm::ESHandle<CaloTPGTranscoder> L1EGCrystalClusterEmulatorProducer::decoder_

private

Definition at line 242 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

ebGeometry

const CaloSubdetectorGeometry* L1EGCrystalClusterEmulatorProducer::ebGeometry

private

Definition at line 247 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

ecalTPEBToken_

edm::EDGetTokenT<EcalEBTrigPrimDigiCollection> L1EGCrystalClusterEmulatorProducer::ecalTPEBToken_

private

Definition at line 240 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

hbGeometry

const CaloSubdetectorGeometry* L1EGCrystalClusterEmulatorProducer::hbGeometry

private

Definition at line 248 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

hbTopology

edm::ESHandle<HcalTopology> L1EGCrystalClusterEmulatorProducer::hbTopology

private

Definition at line 249 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

hcalTPToken_

edm::EDGetTokenT<edm::SortedCollection<HcalTriggerPrimitiveDigi> > L1EGCrystalClusterEmulatorProducer::hcalTPToken_

private

Definition at line 241 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

hcTopology_

const HcalTopology* L1EGCrystalClusterEmulatorProducer::hcTopology_

private

Definition at line 250 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().