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

	EDProducer (const EDProducer &)=delete

bool	hasAbilityToProduceInBeginLumis () const final

bool	hasAbilityToProduceInBeginProcessBlocks () const final

bool	hasAbilityToProduceInBeginRuns () const final

bool	hasAbilityToProduceInEndLumis () const final

bool	hasAbilityToProduceInEndProcessBlocks () const final

bool	hasAbilityToProduceInEndRuns () const final

const EDProducer &	operator= (const EDProducer &)=delete

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

Additional Inherited Members
Public Types inherited from edm::stream::EDProducer<>
using	CacheTypes = CacheContexts< T... >

using	GlobalCache = typename CacheTypes::GlobalCache

using	HasAbility = AbilityChecker< T... >

using	InputProcessBlockCache = typename CacheTypes::InputProcessBlockCache

using	LuminosityBlockCache = typename CacheTypes::LuminosityBlockCache

using	LuminosityBlockContext = LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >

using	LuminosityBlockSummaryCache = typename CacheTypes::LuminosityBlockSummaryCache

using	RunCache = typename CacheTypes::RunCache

using	RunContext = RunContextT< RunCache, GlobalCache >

using	RunSummaryCache = typename CacheTypes::RunSummaryCache

Detailed Description

Definition at line 227 of file L1EGammaCrystalsEmulatorProducer.cc.

Constructor & Destructor Documentation

◆ L1EGCrystalClusterEmulatorProducer()

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

explicit

Definition at line 340 of file L1EGammaCrystalsEmulatorProducer.cc.

     : ecalTPEBToken_(consumes<EcalEBTrigPrimDigiCollection>(iConfig.getParameter<edm::InputTag>("ecalTPEB"))),
       hcalTPToken_(
           consumes<edm::SortedCollection<HcalTriggerPrimitiveDigi> >(iConfig.getParameter<edm::InputTag>("hcalTP"))),
       decoderTag_(esConsumes<CaloTPGTranscoder, CaloTPGRecord>(edm::ESInputTag("", ""))),
       calib_(iConfig.getParameter<edm::ParameterSet>("calib")),
       caloGeometryTag_(esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag("", ""))),
       hbTopologyTag_(esConsumes<HcalTopology, HcalRecNumberingRecord>(edm::ESInputTag("", ""))) {
   produces<l1tp2::CaloCrystalClusterCollection>();
   produces<BXVector<l1t::EGamma> >();
   produces<l1tp2::CaloTowerCollection>("L1CaloTowerCollection");
 }

◆ ~L1EGCrystalClusterEmulatorProducer()

L1EGCrystalClusterEmulatorProducer::~L1EGCrystalClusterEmulatorProducer ( )

override

Definition at line 353 of file L1EGammaCrystalsEmulatorProducer.cc.

353 {}

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 1183 of file L1EGammaCrystalsEmulatorProducer.cc.

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

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

Referenced by produce().

◆ passes_looseTkiso()

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

private

Definition at line 1197 of file L1EGammaCrystalsEmulatorProducer.cc.

                                                                               {
   bool is_iso = (b0 + b1 * std::exp(-b2 * pt) > iso);
   if (pt > plateau_ss)
     is_iso = true;
   return is_iso;
 }

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

Referenced by produce().

◆ passes_looseTkss()

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

private

Definition at line 1218 of file L1EGammaCrystalsEmulatorProducer.cc.

                                                                             {
   bool is_ss = ((e0_looseTkss - e1_looseTkss * std::exp(-e2_looseTkss * pt)) <= ss);
   if (pt > plateau_ss)
     is_ss = true;
   return is_ss;
 }

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

Referenced by produce().

◆ passes_photon()

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

private

Definition at line 1211 of file L1EGammaCrystalsEmulatorProducer.cc.

                                                                           {
   bool is_ss = (pss > d0 - d1 * pt);
   if (pt > plateau_ss)
     is_ss = true;
   return is_ss;
 }

References d0, d1, plateau_ss, and DiDispStaMuonMonitor_cfi::pt.

Referenced by produce().

◆ passes_ss()

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

private

Definition at line 1204 of file L1EGammaCrystalsEmulatorProducer.cc.

                                                                      {
   bool is_ss = ((c0_ss + c1_ss * std::exp(-c2_ss * pt)) <= ss);
   if (pt > plateau_ss)
     is_ss = true;
   return is_ss;
 }

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

Referenced by produce().

◆ produce()

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

overrideprivate

Definition at line 355 of file L1EGammaCrystalsEmulatorProducer.cc.

                                                                                               {
   using namespace edm;
  
   edm::Handle<EcalEBTrigPrimDigiCollection> pcalohits;
   iEvent.getByToken(ecalTPEBToken_, pcalohits);
  
   const auto& caloGeometry = iSetup.getData(caloGeometryTag_);
   ebGeometry = caloGeometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
   hbGeometry = caloGeometry.getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
   const auto& hbTopology = iSetup.getData(hbTopologyTag_);
   hcTopology_ = &hbTopology;
   HcalTrigTowerGeometry theTrigTowerGeometry(hcTopology_);
  
   const auto& decoder = iSetup.getData(decoderTag_);
  
   //****************************************************************
   //******************* 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), "L1CaloTowerCollection");
 }

Member Data Documentation

◆ calib_

l1tp2::ParametricCalibration L1EGCrystalClusterEmulatorProducer::calib_

private

Definition at line 245 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

◆ caloGeometryTag_

edm::ESGetToken<CaloGeometry, CaloGeometryRecord> L1EGCrystalClusterEmulatorProducer::caloGeometryTag_

private

Definition at line 247 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

◆ decoderTag_

edm::ESGetToken<CaloTPGTranscoder, CaloTPGRecord> L1EGCrystalClusterEmulatorProducer::decoderTag_

private

Definition at line 243 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

◆ ebGeometry

const CaloSubdetectorGeometry* L1EGCrystalClusterEmulatorProducer::ebGeometry

private

Definition at line 248 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

◆ ecalTPEBToken_

edm::EDGetTokenT<EcalEBTrigPrimDigiCollection> L1EGCrystalClusterEmulatorProducer::ecalTPEBToken_

private

Definition at line 241 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

◆ hbGeometry

const CaloSubdetectorGeometry* L1EGCrystalClusterEmulatorProducer::hbGeometry

private

Definition at line 249 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

◆ hbTopologyTag_

edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> L1EGCrystalClusterEmulatorProducer::hbTopologyTag_

private

Definition at line 250 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

◆ hcalTPToken_

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

private

Definition at line 242 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

◆ hcTopology_

const HcalTopology* L1EGCrystalClusterEmulatorProducer::hcTopology_

private

Definition at line 251 of file L1EGammaCrystalsEmulatorProducer.cc.

Referenced by produce().

Private Attributes
l1tp2::ParametricCalibration	calib_

edm::ESGetToken< CaloGeometry, CaloGeometryRecord >	caloGeometryTag_

edm::ESGetToken< CaloTPGTranscoder, CaloTPGRecord >	decoderTag_

const CaloSubdetectorGeometry *	ebGeometry

edm::EDGetTokenT< EcalEBTrigPrimDigiCollection >	ecalTPEBToken_

const CaloSubdetectorGeometry *	hbGeometry

edm::ESGetToken< HcalTopology, HcalRecNumberingRecord >	hbTopologyTag_

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

const HcalTopology *	hcTopology_

Classes

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ L1EGCrystalClusterEmulatorProducer()

◆ ~L1EGCrystalClusterEmulatorProducer()

Member Function Documentation

◆ get_calibrate()

◆ passes_iso()

◆ passes_looseTkiso()

◆ passes_looseTkss()

◆ passes_photon()

◆ passes_ss()

◆ produce()

Member Data Documentation

◆ calib_

◆ caloGeometryTag_

◆ decoderTag_

◆ ebGeometry

◆ ecalTPEBToken_

◆ hbGeometry

◆ hbTopologyTag_

◆ hcalTPToken_

◆ hcTopology_