eFexTowerBuilder redesign, and new algorithm to convert eFexTowers to eTowers in the L1CaloFEX simulation

• eFexTowerBuilder redesign that will take the mapping from a given file, or otherwise try to generate it in the first event.
• eTowerMakerFromEfexTowers algorithm that produces a collection of eTowers from eFexTowers

The mapping file for the eFexTowerBuilder needs adding to the calib file area so I should do that before we accept this MR so that I can update the path in the algo property. Link to calibarea upload request: https://its.cern.ch/jira/browse/ATLASG-2463

Trigger/TrigT1/L1CaloFEX/L1CaloFEXAlgos/src/eFexTowerBuilder.cxx

@@ -23,6 +23,7 @@
 
 #include "TFile.h"
 #include "TTree.h"
+#include "PathResolver/PathResolver.h"
 
 
 namespace LVL1 {
@@ -40,57 +41,138 @@ StatusCode eFexTowerBuilder::initialize() {
     CHECK( m_scellKey.initialize(true) );
     CHECK( m_outKey.initialize(true) );
 
+    if (auto fileName = PathResolverFindCalibFile( m_mappingFile ); !fileName.empty()) {
+        std::unique_ptr<TFile> f( TFile::Open(fileName.c_str()) );
+        if (f) {
+            TTree* t = f->Get<TTree>("mapping");
+            if(t) {
+                unsigned long long scid = 0;
+                std::pair<int,int> coord = {0,0};
+                std::pair<int,int> slot;
+                t->SetBranchAddress("scid",&scid);
+                t->SetBranchAddress("etaIndex",&coord.first);
+                t->SetBranchAddress("phiIndex",&coord.second);
+                t->SetBranchAddress("slot1",&slot.first);
+                t->SetBranchAddress("slot2",&slot.second);
+                for(Long64_t i=0;i<t->GetEntries();i++) {
+                    t->GetEntry(i);
+                    m_scMap[scid] = std::make_pair(coord,slot);
+                }
+            }
+        }
+        if (m_scMap.empty()) {
+            ATH_MSG_WARNING("Failed to load sc -> eFexTower map from " << fileName);
+        } else {
+            ATH_MSG_INFO("Loaded sc -> eFexTower map from " << fileName);
+        }
+    }
 
 
     return StatusCode::SUCCESS;
 }
 
+StatusCode eFexTowerBuilder::fillTowers(const EventContext& ctx) const {
 
-StatusCode eFexTowerBuilder::execute(const EventContext& ctx) const {
-    ATH_MSG_DEBUG("Executing " << name() << "...");
-    setFilterPassed(true, ctx);
 
     SG::ReadCondHandle<CaloSuperCellDetDescrManager> ddm{m_ddmKey,ctx};
-
     SG::ReadHandle<xAOD::TriggerTowerContainer> tTowers(m_ttKey,ctx);
     if(!tTowers.isValid()){
         ATH_MSG_FATAL("Could not retrieve collection " << m_ttKey.key() );
         return StatusCode::FAILURE;
     }
-    SG::ReadHandle<CaloCellContainer> scells(m_scellKey,ctx);
+    SG::ReadHandle<CaloCellContainer> scells(m_scellKey,ctx); // n.b. 34048 is a full complement of scells
     if(!scells.isValid()){
         ATH_MSG_FATAL("Could not retrieve collection " << m_scellKey.key() );
         return StatusCode::FAILURE;
     }
 
+    std::map<std::pair<int,int>,std::array<int,11>> towers;
+
+    for (auto digi: *scells) {
+        const auto itr = m_scMap.find(digi->ID().get_compact());
+        if (itr == m_scMap.end()) { continue; } // not in map so not mapping to a tower
+        int val =  std::round(digi->energy()/(12.5*std::cosh(digi->eta())));
+        bool isMasked = ((digi)->provenance()&0x80);
+        if (isMasked) val = std::numeric_limits<int>::max();
+        auto& tower = towers[itr->second.first];
+        if (itr->second.second.second<11) {
+            // doing an energy split between slots ... don't include a masked channel
+            if (!isMasked) {
+                tower.at(itr->second.second.first) += val >> 1;
+                tower.at(itr->second.second.second) += (val - (val >> 1));
+            }
+        } else {
+            tower.at(itr->second.second.first) += val;
+        }
+    }
+
+    // add tile energies from TriggerTowers
+    static const auto etaIndex = [](float eta) { return int( eta*10 ) + ((eta<0) ? -1 : 1); };
+    static const auto phiIndex = [](float phi) { return int( phi*32./M_PI ) + (phi<0 ? -1 : 1); };
+    for(const auto& tTower : *tTowers) {
+        if (std::abs(tTower->eta()) > 1.5) continue;
+        if (tTower->sampling() != 1) continue;
+        double phi = tTower->phi(); if(phi > M_PI) phi -= 2.*M_PI;
+        towers[std::pair(etaIndex(tTower->eta()),phiIndex(phi))][10] = tTower->cpET();
+    }
+
     SG::WriteHandle<xAOD::eFexTowerContainer> eTowers = SG::WriteHandle<xAOD::eFexTowerContainer>(m_outKey,ctx);
     ATH_CHECK( eTowers.record(std::make_unique<xAOD::eFexTowerContainer>(),std::make_unique<xAOD::eFexTowerAuxContainer>()) );
 
-    struct TowerSCells {
-        std::vector<int> ps;
-        std::vector<std::pair<float,int>> l1;
-        std::vector<std::pair<float,int>> l2;
-        std::vector<int> l3;
-        std::vector<int> had;
-        std::vector<int> other;
-    };
-
-    auto calToFex = [&](int calEt) {
-        if (m_mappingVerificationMode.value()) return calEt;
-        if(calEt == std::numeric_limits<int>::max()) return 1024; // indicates masked channel
+    static const auto calToFex = [](int calEt) {
+        if(calEt == std::numeric_limits<int>::max()) return 0; // indicates masked channel
         if(calEt<448) return std::max((calEt&~1)/2+32,1);
         if(calEt<1472) return (calEt-448)/4+256;
         if(calEt<3520) return (calEt-1472)/8+512;
         if(calEt<11584) return (calEt-3520)/32+768;
         return 1020;
     };
-    auto etaIndex = [](float eta) { return int( eta*10 ) + ((eta<0) ? -1 : 1); };
-    auto phiIndex = [](float phi) { return int( phi*32./ROOT::Math::Pi() ) + (phi<0 ? -1 : 1); };
+
+    // now create the towers
+    for(auto& [coord,counts] : towers) {
+        size_t ni = (std::abs(coord.first)<=15) ? 10 : 11; // ensures we skip the tile towers for next line
+        for(size_t i=0;i<ni;++i) counts[i] = calToFex(counts[i]); // do latome energy scaling to non-tile towers
+        eTowers->push_back( std::make_unique<xAOD::eFexTower>() );
+        eTowers->back()->initialize( ( (coord.first<0 ? 0.5:-0.5) + coord.first)*0.1 ,
+                                 ( (coord.second<0 ? 0.5:-0.5) + coord.second)*M_PI/32,
+                                     std::vector<uint16_t>(counts.begin(), counts.end()),
+                                 -1, /* module number */
+                                 -1, /* fpga number */
+                                 0,0 /* status flags ... could use to indicate which cells were actually present?? */);
+
+    }
+
+    return StatusCode::SUCCESS;
+
+}
+
+StatusCode eFexTowerBuilder::fillMap(const EventContext& ctx) const {
+
+    ATH_MSG_INFO("Filling sc -> eFexTower map");
+
+    SG::ReadCondHandle<CaloSuperCellDetDescrManager> ddm{m_ddmKey,ctx};
+    SG::ReadHandle<CaloCellContainer> scells(m_scellKey,ctx); // 34048 is a full complement of scells
+    if(!scells.isValid()){
+        ATH_MSG_FATAL("Could not retrieve collection " << m_scellKey.key() );
+        return StatusCode::FAILURE;
+    }
+    if (scells->size() != 34048) {
+        ATH_MSG_FATAL("Cannot fill sc -> eFexTower mapping with an incomplete sc collection");
+        return StatusCode::FAILURE;
+    }
+    struct TowerSCells {
+        std::vector<unsigned long long> ps;
+        std::vector<std::pair<float,unsigned long long>> l1;
+        std::vector<std::pair<float,unsigned long long>> l2;
+        std::vector<unsigned long long> l3;
+        std::vector<unsigned long long> had;
+        std::vector<unsigned long long> other;
+    };
+    static const auto etaIndex = [](float eta) { return int( eta*10 ) + ((eta<0) ? -1 : 1); };
+    static const auto phiIndex = [](float phi) { return int( phi*32./ROOT::Math::Pi() ) + (phi<0 ? -1 : 1); };
     std::map<std::pair<int,int>,TowerSCells> towers;
-    std::map<std::pair<int,int>,std::set<int>> towerIDs;
-    // sort cells by eta, phi, into layers
-    // deliberately leaving next line commented to restore for debugging in future
-    //TFile f1("test.root","RECREATE");TTree* badChans = new TTree("badChans","badChans");uint32_t offId; badChans->Branch("offId",&offId);uint32_t status; badChans->Branch("status",&status);
+    std::map<unsigned long long,int> eTowerSlots; // not used by this alg, but we produce the map for benefit of eFexTower->eTower alg
+
     for (auto digi: *scells) {
         Identifier id = digi->ID(); // this is if using supercells
 
@@ -98,25 +180,13 @@ StatusCode eFexTowerBuilder::execute(const EventContext& ctx) const {
             float eta = elem->eta_raw(); // this seems more symmetric
             int sampling = elem->getSampling();
             if(sampling==6 && ddm->getCaloCell_ID()->region(id)==0 && eta<0) eta-=0.01; // nudge this L2 endcap supercell into correct tower (right on boundary)
-            auto val =  round(digi->energy()/(12.5*std::cosh(digi->eta())));
 
-            bool isMasked = ((digi)->provenance()&0x80);
+            unsigned long long val = id.get_compact();
+
+            int towerid = -1;int slot = -1;bool issplit = false;
+            CHECK(m_eFEXSuperCellTowerIdProviderTool->geteTowerIDandslot(id.get_compact(), towerid, slot, issplit));
+            eTowerSlots[id.get_compact()] = slot;
 
-            if( isMasked ) {
-                // deliberately leaving following lines commented in case we want to restore for tests
-                //badChans.SetPoint(badChans.GetN(),eta,elem->phi_raw());
-                //offId =  id.get_identifier32().get_compact();
-                //status = bcCont->offlineStatus(id).packedData();
-                //badChans->Fill();
-                val = std::numeric_limits<int>::max();
-            }
-            if(m_mappingVerificationMode.value()) {
-                int towerid = -1;int slot = -1;bool issplit = false;
-                CHECK(m_eFEXSuperCellTowerIdProviderTool->geteTowerIDandslot(id.get_compact(), towerid, slot, issplit));
-                towerIDs[std::pair(etaIndex(eta), phiIndex(elem->phi_raw()))].insert(towerid);
-                // will use slot+1 as the energy value in this mode
-                val = slot+1;
-            }
             auto& sc = towers[std::pair(etaIndex(eta),phiIndex(elem->phi_raw()))];
             switch(sampling) {
                 case 0: case 4: //lar barrel/endcap presampler
@@ -135,16 +205,11 @@ StatusCode eFexTowerBuilder::execute(const EventContext& ctx) const {
             }
         }
     }
-    //badChans->Write("badchans");f1.Close(); // comment left here so we know how bad channels were verified
-    // add tile energies from TriggerTowers
-    for(const auto& tTower : *tTowers) {
-        if (std::abs(tTower->eta()) > 1.5) continue;
-        if (tTower->sampling() != 1) continue;
-        double phi = tTower->phi(); if(phi > ROOT::Math::Pi()) phi -= 2.*ROOT::Math::Pi();
-        towers[std::pair(etaIndex(tTower->eta()),phiIndex(phi))].had.push_back(m_mappingVerificationMode.value() ? 11 : tTower->cpET());
-    }
+
 
     // sort (by increasing eta) l1/l2 sc and handle special cases
+    // finally also output the eTower slot vector
+    std::vector<size_t> slotVector(11);
     for(auto& [coord,sc] : towers) {
         std::sort(sc.l1.begin(),sc.l1.end());
         std::sort(sc.l2.begin(),sc.l2.end());
@@ -166,56 +231,77 @@ StatusCode eFexTowerBuilder::execute(const EventContext& ctx) const {
         }
         // handle case @ |eta|~1.8-2 with 6 L1 cells
         if (sc.l1.size()==6) {
-            sc.l1.at(0).second += int(sc.l1.at(1).second>>1);
-            sc.l1.at(1).second = int(sc.l1.at(1).second>>1) + sc.l1.at(2).second;
-            sc.l1.at(2).second = sc.l1.at(3).second + int(sc.l1.at(4).second>>1);
-            sc.l1.at(3).second = int(sc.l1.at(4).second>>1) + sc.l1.at(5).second;
-            sc.l1.resize(4);
+            m_scMap[sc.l1.at(0).second] = std::pair(coord,std::pair(1,11));
+            m_scMap[sc.l1.at(1).second] = std::pair(coord,std::pair(1,2));
+            m_scMap[sc.l1.at(2).second] = std::pair(coord,std::pair(2,11));
+            m_scMap[sc.l1.at(3).second] = std::pair(coord,std::pair(3,11));
+            m_scMap[sc.l1.at(4).second] = std::pair(coord,std::pair(3,4));
+            m_scMap[sc.l1.at(5).second] = std::pair(coord,std::pair(4,11));
+            slotVector[1] = eTowerSlots[sc.l1.at(0).second];
+            slotVector[2] = eTowerSlots[sc.l1.at(2).second];
+            slotVector[3] = eTowerSlots[sc.l1.at(3).second];
+            slotVector[4] = eTowerSlots[sc.l1.at(5).second];
         }
 
         // for |eta|>2.4 there's only 1 l1 sc, to match hardware this should be compared placed in the 'last' l1 input
         if (sc.l1.size()==1) {
-            auto tmp = sc.l1.at(0); sc.l1.clear(); sc.l1.resize(4,std::pair(0,0)); sc.l1.at(3) = tmp;
+            m_scMap[sc.l1.at(0).second] = std::pair(coord,std::pair(4,11));
+            slotVector[1] = 1; slotVector[2] = 2; slotVector[3] = 3; slotVector[4] = eTowerSlots[sc.l1.at(0).second];
         }
 
-        // create an eFexTower ... use module number = -1, fpga = -1
-        std::vector<uint16_t> counts; counts.reserve(11);
-        counts.push_back(sc.ps.empty() ? 0 : calToFex(sc.ps.at(0)));
-        for(size_t i=0;i<4;i++) counts.push_back((sc.l1.size() > i) ? calToFex(sc.l1.at(i).second) : 0);
-        for(size_t i=0;i<4;i++) counts.push_back((sc.l2.size() > i) ? calToFex(sc.l2.at(i).second) : 0);
-        counts.push_back(sc.l3.empty() ? 0 : calToFex(sc.l3.at(0)));
-        counts.push_back(sc.had.empty() ? 0 : (std::abs(coord.first) <= 15 ? sc.had.at(0) : calToFex(sc.had.at(0))) ); // tile needs no convert
-
-
-        if(m_mappingVerificationMode.value()) {
-            std::cout << coord.first << " " << coord.second << " : ";
-            for(auto& id : towerIDs[coord]) std::cout << id << ",";
-            std::cout << " : ";
-            // counts(-1) should typically be just 0,1,2,3,...
-            // -1 indicates no input
-            // exceptions are:
-            //   |eta|=25 : the 4th l1 will be a "0" (treating as ps)
-            //   |eta|=18,19,20 :
-            for(size_t i=0;i<counts.size();i++) {
-                std::cout << counts[i]-1 << ",";
-            }
-            std::cout << std::endl;
-            continue;
-        }
-
-        eTowers->push_back( std::make_unique<xAOD::eFexTower>() );
-        eTowers->back()->initialize( ( (coord.first<0 ? 0.5:-0.5) + coord.first)*0.1 ,
-                                     ( (coord.second<0 ? 0.5:-0.5) + coord.second)*ROOT::Math::Pi()/32,
-                                     counts,
-                                     -1, /* module number */
-                                     -1, /* fpga number */
-                                     0,0 /* status flags ... could use to indicate which cells were actually present?? */);
-
+        // fill the map with sc ids -> tower coord + slot
+        if (!sc.ps.empty()) {m_scMap[sc.ps.at(0)] = std::pair(coord,std::pair(0,11)); slotVector[0] = eTowerSlots[sc.ps.at(0)]; }
+        if(sc.l1.size()==4) for(size_t i=0;i<4;i++) if(sc.l1.size() > i) {m_scMap[sc.l1.at(i).second] = std::pair(coord,std::pair(i+1,11)); slotVector[i+1] = eTowerSlots[sc.l1.at(i).second]; }
+        for(size_t i=0;i<4;i++) if(sc.l2.size() > i) { m_scMap[sc.l2.at(i).second] = std::pair(coord,std::pair(i+5,11)); slotVector[i+5] = eTowerSlots[sc.l2.at(i).second]; }
+        if (!sc.l3.empty()) {m_scMap[sc.l3.at(0)] = std::pair(coord,std::pair(9,11)); slotVector[9] = eTowerSlots[sc.l3.at(0)]; }
+        if (!sc.had.empty()) {m_scMap[sc.had.at(0)] = std::pair(coord,std::pair(10,11));slotVector[10] = eTowerSlots[sc.had.at(0)]; }
+
+        // finally output the slotVector for this tower
+        // do only for the slots that don't match
+        // note to self: seems like everything is fine apart from the l1->ps remap for |eta|>2.4
+        // so leaving this bit commented out for now ... useful to leave it here in case need to recheck in future
+//        for(size_t i=0;i<slotVector.size();i++) {
+//            if(slotVector[i] != i) {
+//                std::cout << coord.first << "," << coord.second << "," << i << "," << slotVector[i] << std::endl;
+//            }
+//        }
     }
 
+    // save the map to disk
+    TFile f("scToEfexTowers.root","RECREATE");
+    TTree* t = new TTree("mapping","mapping");
+    unsigned long long scid = 0;
+    std::pair<int,int> coord = {0,0};
+    std::pair<int,int> slot = {-1,-1};
+    t->Branch("scid",&scid);
+    t->Branch("etaIndex",&coord.first);
+    t->Branch("phiIndex",&coord.second);
+    t->Branch("slot1",&slot.first);
+    t->Branch("slot2",&slot.second);
+    for(auto& [id,val] : m_scMap) {
+        scid = id; coord = val.first; slot = val.second;
+        t->Fill();
+    }
+    t->Write();
+    f.Close();
 
     return StatusCode::SUCCESS;
 
 }
 
+
+StatusCode eFexTowerBuilder::execute(const EventContext& ctx) const {
+    ATH_MSG_DEBUG("Executing " << name() << "...");
+    setFilterPassed(true, ctx);
+
+
+    {
+        std::lock_guard lock(m_fillMapMutex);
+        if (m_scMap.empty()) CHECK( fillMap(ctx) );
+    }
+
+    return fillTowers(ctx);
+
+}
+
 } // LVL1 Namespace
\ No newline at end of file