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	hasAbilityToProduceInBeginRuns () const final
bool	hasAbilityToProduceInEndLumis () 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, begin, 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_, 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().