From 4bad9efc9c135a30fbd90e9ce5b3a8b3c2965742 Mon Sep 17 00:00:00 2001 From: amyers Date: Mon, 17 Jun 2019 18:03:59 +0200 Subject: [PATCH 1/5] Migrating changes from l1calo-run3-offline --- .../TrigT1CaloFexSim/JGTowerReader.h | 2 + .../TrigT1CaloFexSim/TrigT1CaloFexSim/JwoJ.h | 64 ++- .../TrigT1CaloFexSim/METAlg.h | 15 +- .../TrigT1CaloFexSim/TrigT1CaloFexSim/Rho.h | 123 ++---- .../TrigT1CaloFexSim/src/JGTowerReader.cxx | 379 +++++++++++++++++- .../TrigT1/TrigT1CaloFexSim/src/METAlg.cxx | 332 ++++++++++----- 6 files changed, 687 insertions(+), 228 deletions(-) diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JGTowerReader.h b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JGTowerReader.h index db4343368ca..e2ccd7d12a7 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JGTowerReader.h +++ b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JGTowerReader.h @@ -30,12 +30,14 @@ #include "TH2.h" #include "TrigT1CaloFexSim/JetAlg.h" #include "TrigT1CaloFexSim/METAlg.h" + #include "string.h" std::vector splitString(std::string parentString, std::string sep, bool stripEmpty=false); class JGTowerReader: public ::AthAlgorithm { public: JGTowerReader( const std::string& name, ISvcLocator* pSvcLocator ); virtual ~JGTowerReader(); + virtual StatusCode MET_etaBins(const xAOD::JGTowerContainer* gTs, TString metName, bool usegFEX, bool useRhoSub, bool useJwoJ, bool usePUfit); virtual StatusCode initialize(); virtual StatusCode execute(); virtual StatusCode finalize(); diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JwoJ.h b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JwoJ.h index bb96f74ccef..54be2b5e8fc 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JwoJ.h +++ b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JwoJ.h @@ -35,6 +35,8 @@ void BuildBlocksFromTowers(std::vector& blocks, const xAOD:: std::vector neighbors = grid.neighbors(*seed, blockRows, blockCols); float seed_Et = seed->et(); + float eta = fabs(seed->eta()); + if(!useNegTowers) seed_Et = TMath::Abs(seed_Et); double block_area(0.0); double block_pt(seed_Et); @@ -44,6 +46,8 @@ void BuildBlocksFromTowers(std::vector& blocks, const xAOD:: const xAOD::JGTower* neighbor = towers.at(neighborIndex); block_area += neighbor->deta()*neighbor->dphi(); neighbor_pt = neighbor->et(); + float n_eta = fabs(neighbor->eta()); + if(!useNegTowers) neighbor_pt = TMath::Abs(neighbor_pt); block_pt += neighbor_pt; } @@ -65,7 +69,7 @@ void BuildBlocksFromTowers(std::vector& blocks, const xAOD:: *@brief Separates Et into hard and soft terms depending on threshold pt. Stores hard and soft MET terms in a vector to be passed. *@return @c std::vector */ -std::vector Run_JwoJ(const xAOD::JGTowerContainer* towers, float pTcone_cut, bool useNegTowers){ +std::vector Run_JwoJ(const xAOD::JGTowerContainer* towers, float rho, float pTcone_cut, bool useEtaBins, bool useNegTowers){ std::vector blocks; BuildBlocksFromTowers(blocks, *towers, 3, 3, useNegTowers); @@ -83,10 +87,14 @@ std::vector Run_JwoJ(const xAOD::JGTowerContainer* towers, float pTcone_c float pt_cone = blocks[b].Pt(); float seed_Et = (towers->at(blocks[b].seedIndex()))->et(); - float block_etx = seed_Et*cos(block_phi); - float block_ety = seed_Et*sin(block_phi); + float block_etx = 0; + float block_ety = 0; - if(TMath::Abs(blocks[b].Eta()) < 2.4){ + if(useEtaBins){ + seed_Et -= rho; + block_etx = seed_Et*TMath::Cos(block_phi); + block_ety = seed_Et*TMath::Sin(block_phi); + if(pt_cone > pTcone_cut){ Ht += seed_Et; Htx += block_etx; @@ -95,17 +103,37 @@ std::vector Run_JwoJ(const xAOD::JGTowerContainer* towers, float pTcone_c Et += seed_Et; Etx += block_etx; Ety += block_ety; - } - }else{ - if(pt_cone > 4*pTcone_cut){ - Ht += seed_Et; - Htx += block_etx; - Hty += block_ety; } - else{ - Et += seed_Et; - Etx += block_etx; - Ety += block_ety; + }else{ + if(TMath::Abs(blocks[b].Eta()) < 2.4){ + seed_Et -= rho; + block_etx = seed_Et*TMath::Cos(block_phi); + block_ety = seed_Et*TMath::Sin(block_phi); + + if(pt_cone > pTcone_cut){ + Ht += seed_Et; + Htx += block_etx; + Hty += block_ety; + }else{ + Et += seed_Et; + Etx += block_etx; + Ety += block_ety; + } + }else{ + seed_Et -= 4*rho; + block_etx = seed_Et*TMath::Cos(block_phi); + block_ety = seed_Et*TMath::Sin(block_phi); + + if(pt_cone > 4*pTcone_cut){ + Ht += seed_Et; + Htx += block_etx; + Hty += block_ety; + } + else{ + Et += seed_Et; + Etx += block_etx; + Ety += block_ety; + } } } scalar_Et += seed_Et; @@ -118,10 +146,14 @@ std::vector Run_JwoJ(const xAOD::JGTowerContainer* towers, float pTcone_c //convert to GeV for fitting process Et_vals.push_back(scalar_Et); - Et_vals.push_back(MHT); + Et_vals.push_back(Htx); + Et_vals.push_back(Hty); + Et_vals.push_back(Etx); + Et_vals.push_back(Ety); + Et_vals.push_back(MHT); Et_vals.push_back(MET); Et_vals.push_back(Et_tot); - + return Et_vals; } diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/METAlg.h b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/METAlg.h index f1f2ecca895..15bce7c7ff2 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/METAlg.h +++ b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/METAlg.h @@ -42,6 +42,11 @@ class METAlg{ float rho = 0; float mht = 0; float mst = 0; + float mht_x = 0; + float mht_y = 0; + float mst_x = 0; + float mst_y = 0; + float scalar_Et = 0; }; static std::map> m_METMap; @@ -52,7 +57,7 @@ class METAlg{ /** *@brief Calculates MET with pileup subtraction */ - static StatusCode SubtractRho_MET(const xAOD::JGTowerContainer* towers, TString metname, bool useRMS, bool useMedian, bool useNegTowers); + static StatusCode SubtractRho_MET(const xAOD::JGTowerContainer* towers, TString metname, bool useEtaBins, bool useRMS, bool useMedian, bool useNegTowers); /** *@brief Calculates MET with Softkiller */ @@ -60,16 +65,14 @@ class METAlg{ /** *@brief Calculates MET with Jets without Jets */ - static StatusCode JwoJ_MET(const xAOD::JGTowerContainer* towers, TString metname, float pTcone_cut, bool useNegTowers); + static StatusCode JwoJ_MET(const xAOD::JGTowerContainer* towers, TString metname, float pTcone_cut, bool useEtaBins, bool useRho, bool useNegTowers); /** *@brief Calculates MET using PUfit */ static StatusCode Pufit_MET(const xAOD::JGTowerContainer* towers, TString metname, bool useNegTowers); - /** - *@brief Calculates MET in bins of eta(and phi) depending on jFEX/gFEX geometry - */ - static StatusCode MET_etaBins(const xAOD::JGTowerContainer* towers, TString metName,bool usegFEX, bool useRhoSub, bool usePUfit); + static float Rho_avg(const xAOD::JGTowerContainer* towers, bool useNegTowers); + static float Rho_avg_etaRings(const xAOD::JGTowerContainer* towers, int fpga, bool useNegTowers); }; #endif diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/Rho.h b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/Rho.h index a8a07d90687..5e46552ce77 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/Rho.h +++ b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/Rho.h @@ -41,60 +41,63 @@ Rho_med(): Calculates rho as the median tower energy density in the barrel, agai *@brief Calculates rho as the average tower energy density in the barrel region of the calorimeter, to be scaled with are in more forward regions. An upper threshold of 10 GeV is applied to eliminate bias from hard scatter events *@return @c float */ -float Rho_bar(const xAOD::JGTowerContainer* towers, const bool useNegTowers){ +float Rho_bar(const xAOD::JGTowerContainer* towers, const bool useEtaBins, int fpga, const bool useNegTowers){ + //fpga = 0 ->entire barrel + //fpga = 1 -> FPGA_A + //fpga = 2 -> FPGA_B + //fpga = 3 -> FPGA_C float rho = 0; float et_max = 10*Gaudi::Units::GeV; //an upper threshold such that the average rho is not biased by hard scatter events const unsigned int size = towers->size(); - int length = 0; + int length = 0; + float area_a = 0; + float area_b = 0; + float area_c = 0; for(unsigned int i = 0; i < size; i++){ const xAOD::JGTower* tower = towers->at(i); + float eta = TMath::Abs(tower->eta()); float Et = tower->et(); if(!useNegTowers && Et < 0) continue; - if(eta < 2.4){ + if(eta < 0 && eta > -2.5){ //FPGA A bounds + if(eta > -2.4) area_a += 1; + else area_a += 0.5; + } + if(eta < 2.5 && eta > 0){ + if(eta < 2.4) area_b += 1.; + else area_b += 0.5; + } + if(fabs(eta) > 2.5){ + if(eta < 3.2) area_c += 1.; + else area_c += 4.; + } + + if(useEtaBins){ if(Et < et_max){ rho += Et; - length++; + //length++; + } + }else{ + if(eta < 2.4){ + if(Et < et_max){ + rho += Et; + length++; + } } } } float rho_bar = rho/length; + if(fpga == 1) rho_bar = rho/area_a; + if(fpga == 2) rho_bar = rho/area_b; + if(fpga ==3) rho_bar = rho/area_c; return rho_bar; } -float Rho_bar(std::vector towers, const bool useNegTowers){ - - float rho = 0; - float et_max = 10*Gaudi::Units::GeV; //an upper threshold such that the average rho is not biased by hard scatter events - - const unsigned int size = towers.size(); - int length = 0; - - for(unsigned int i = 0; i < size; i++){ - const xAOD::JGTower* tower = towers[i]; - float eta = TMath::Abs(tower->eta()); - float Et = tower->et(); - - if(!useNegTowers && Et < 0) continue; - - if(eta < 2.4){ - if(Et < et_max){ - rho += Et; - length++; - } - } - } - float rho_bar = length > 0 ? rho/length : 0; - - return rho_bar; -} - - /** *@brief Calculates rho as the median tower energy density in the barrel region of the calorimeter, to be scaled with area in more forward regions *@return @c float @@ -123,60 +126,4 @@ float Rho_med(const xAOD::JGTowerContainer* towers, const bool useNegTowers){ return rho; } -const std::vector rhoSub(std::vector towers, bool useNegTowers){ - float EtMiss = 0; - float Ex = 0, Ey = 0, Ex_ = 0, Ey_ = 0; - float threshold = 0; - std::vector met; - - //can calculate rho as either the average or median gTower energy in the barrel - float rho = Rho_bar(towers, false); - - unsigned int size = towers.size(); - - TH1F* h_Et = new TH1F("h_Et", "", 50, 0, 5000); - - for(unsigned int t = 0; t < size; t++){ - const xAOD::JGTower* tower = towers[t]; - float Et = tower->et(); - if(!useNegTowers && Et < 0) continue; - h_Et->Fill(Et); - } - threshold = 3*h_Et->GetRMS(); - - for(unsigned int t = 0; t < size; t++){ - const xAOD::JGTower* tower = towers[t]; - float Et = tower->et(); - float phi = tower->phi(); - float eta = TMath::Abs(tower->eta()); - - float Et_sub = 0; - if(eta < 2.4) Et_sub = TMath::Abs(Et) - rho; - else Et_sub = TMath::Abs(Et) - 4*rho; - - if(Et_sub < threshold) continue; - - //if tower originally had negative Et, keep Et_sub negative - if(useNegTowers && Et < 0) Et_sub = -Et_sub; - - Ex_ = Et_sub*TMath::Cos(phi); - Ey_ = Et_sub*TMath::Sin(phi); - - Ex += Ex_; - Ey += Ey_; - } - - EtMiss = TMath::Sqrt(Ex*Ex + Ey*Ey); - float phi_met = 0; - if(EtMiss != 0) phi_met = TMath::ACos(Ex/EtMiss); - if (Ey<0) phi_met = -phi_met; - - met.push_back(EtMiss); - met.push_back(phi_met); - - delete h_Et; - - return met; -} - #endif diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx b/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx index 6655fd6740f..da2bbd8cac0 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx +++ b/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx @@ -1,5 +1,4 @@ -/* - * Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +/* * Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration */ // TrigT1CaloFexSim includes @@ -33,6 +32,11 @@ #include "xAODTrigger/EnergySumRoIAuxInfo.h" #include "TFile.h" #include "PathResolver/PathResolver.h" +#include "AthContainers/DataVector.h" +#include "AthContainers/ConstDataVector.h" + +#include "CaloTriggerTool/GTowerSCMap.h" +#include "CaloTriggerTool/JTowerSCMap.h" JGTowerReader::JGTowerReader( const std::string& name, ISvcLocator* pSvcLocator ) : @@ -110,6 +114,8 @@ JGTowerReader::~JGTowerReader() { delete jJetSeeds; delete gSeeds; delete acc_rho; + delete acc_mht; + delete acc_mst; METAlg::m_METMap.clear(); jL1Jets.clear(); @@ -149,6 +155,7 @@ StatusCode JGTowerReader::initialize() { acc_mht = new SG::AuxElement::Accessor("MHT"); acc_mst = new SG::AuxElement::Accessor("MST"); + return StatusCode::SUCCESS; } @@ -160,7 +167,15 @@ StatusCode JGTowerReader::finalize() { StatusCode JGTowerReader::execute() { - ATH_MSG_DEBUG ("Executing " << name() << " on event " << m_eventCount); + if(m_eventCount % 100 == 0){ + ATH_MSG_INFO ("Executing " << name() << " on event " << m_eventCount); + } + else if(m_eventCount % 10 == 0){ + ATH_MSG_INFO (" Executing " << name() << " on event " << m_eventCount); + } + else { + ATH_MSG_DEBUG ("Executing " << name() << " on event " << m_eventCount); + } m_eventCount += 1; setFilterPassed(false); //optional: start with algorithm not passed @@ -189,11 +204,107 @@ StatusCode JGTowerReader::execute() { const xAOD::JGTowerContainer* gTowers =0; CHECK( evtStore()->retrieve( gTowers,"GTower")); + //combine EM and Had layers + //ConstDataVector gT_combined(SG::VIEW_ELEMENTS); + + //std::cout<<"Collapsing EM and Had layers"<setStore(gCaloTowersAux); + for(unsigned t = 0; t < gTowers->size(); t++){ + const xAOD::JGTower* gt_em = gTowers->at(t); + const float eta = gt_em->eta(); + const float phi = gt_em->phi(); + const float deta = gt_em->deta(); + const float dphi = gt_em->dphi(); + float totalEt = 0; + + if(fabs(eta) < 3.15 && gt_em->sampling() == 0){ + const xAOD::JGTower* gt_had = gTowers->at(t+544); + totalEt = gt_em->et() + gt_had->et(); + } + //if(t > 544) totalEt = 0; + const float tEt = totalEt; + + const std::vector index(2, 0); + + xAOD::JGTower* m_tower = new xAOD::JGTower(); + gCaloTowers->push_back(m_tower); + const int t_ = t; + m_tower->initialize(t_, eta, phi); + m_tower->setdEta(deta); + m_tower->setdPhi(dphi); + m_tower->setEt(tEt); + m_tower->setSCIndex(index); + m_tower->setTileIndex(index); + const int sampling = gt_em->sampling(); + m_tower->setSampling(sampling); + + } + ConstDataVector temp_a(SG::VIEW_ELEMENTS); + ConstDataVector temp_b(SG::VIEW_ELEMENTS); + ConstDataVector temp_c(SG::VIEW_ELEMENTS); + + for(unsigned int t = 0; t < gCaloTowers->size(); t++){ + const xAOD::JGTower* tower = gCaloTowers->at(t); + float eta = tower->eta(); + + if(eta > -2.5 && eta < 0) temp_a.push_back(tower); + if(eta > 0 && eta < 2.5) temp_b.push_back(tower); + if(fabs(eta) >= 2.5) temp_c.push_back(tower); + } + const xAOD::JGTowerContainer* fpga_a = temp_a.asDataVector(); + const xAOD::JGTowerContainer* fpga_b = temp_b.asDataVector(); + const xAOD::JGTowerContainer* fpga_c = temp_c.asDataVector(); + + xAOD::JGTowerAuxContainer* pu_subAux = new xAOD::JGTowerAuxContainer(); + xAOD::JGTowerContainer* pu_sub = new xAOD::JGTowerContainer(); + pu_sub->setStore(pu_subAux); + + float rhoA = METAlg::Rho_avg_etaRings(fpga_a, 1, false); + float rhoB = METAlg::Rho_avg_etaRings(fpga_b, 2, false); + float rhoC = METAlg::Rho_avg_etaRings(fpga_c, 3, false); + + for(unsigned int t = 0; t < gCaloTowers->size(); t++){ + const xAOD::JGTower* tower = gCaloTowers->at(t); + const float eta = tower->eta(); + const float phi = tower->phi(); + const float deta = tower->deta(); + const float dphi = tower->dphi(); + float Et_sub = tower->et(); + + if(eta < 0 && eta >= -2.5){ + if(eta > -2.4) Et_sub -= rhoA; + else Et_sub -= 0.5*rhoA; + } + if(eta > 0 && eta <= 2.5){ + if(eta < 2.4) Et_sub -= rhoB; + else Et_sub -= 0.5*rhoB; + } + if(fabs(eta) > 2.5){ + if(fabs(eta) < 3.2) Et_sub -= rhoC; + else Et_sub -= 4*rhoC; + } + + const std::vector index(2,0); + const int sampling = tower->sampling(); + xAOD::JGTower* m_tower = new xAOD::JGTower(); + pu_sub->push_back(m_tower); + m_tower->initialize(t, eta, phi); + m_tower->setdEta(deta); + m_tower->setdPhi(dphi); + m_tower->setEt(Et_sub); + m_tower->setSCIndex(index); + m_tower->setTileIndex(index); + m_tower->setSampling(sampling); + } + + //const xAOD::JGTowerContainer* gTs = gT_combined.asDataVector(); //when noise file is not available, set the noise as constant for EM, Had, and FCal respectively if(gT_noise.size()==0){ ATH_MSG_INFO (" Failed to find noise file, setting gTower noise manually "); - for(unsigned int i = 0; i < gTowers->size(); i++) { - gT_noise.push_back(1500); + for(unsigned int i = 0; i < gTowers->size(); i++) { + gT_noise.push_back(1500); } } @@ -211,7 +322,8 @@ StatusCode JGTowerReader::execute() { ATH_MSG_DEBUG ("About to call JFexAlg"); CHECK(JFexAlg(jTowers)); // all the functions for JFex shall be executed here ATH_MSG_DEBUG ("About to call GFexAlg"); - CHECK(GFexAlg(gTowers)); // all the functions for GFex shall be executed here + CHECK(GFexAlg(gCaloTowers)); // all the functions for GFex shall be executed here + //CHECK(GFexAlg(pu_sub)); //run all gFEX algs with pileup subtracted towers ATH_MSG_DEBUG ("About to call ProcessObject()"); CHECK(ProcessObjects()); // this is the function to make output as well as memory cleaning @@ -335,29 +447,264 @@ StatusCode JGTowerReader::GFexAlg(const xAOD::JGTowerContainer* gTs){ CHECK(METAlg::Baseline_MET(gTs,Form("gBaselineNeg%d",NegTowers),gT_noise, NegTowers)); //basic MET reconstruction with a 4 sigma noise cut applied //CHECK(METAlg::SubtractRho_MET(gTs,Form("RhoSubNeg%d",NegTowers),m_useRMS,m_useMedian,NegTowers)); //pileup subtracted MET, can apply dynamic noise cut and use either median or avg rho //No RMS - CHECK(METAlg::SubtractRho_MET(gTs,Form("RhoSubNeg%dMedian",NegTowers),0,1,NegTowers));//Median - CHECK(METAlg::SubtractRho_MET(gTs,Form("RhoSubNeg%dMean",NegTowers),0,0,NegTowers));//Mean + //CHECK(METAlg::SubtractRho_MET(gTs,Form("RhoSubNeg%dMedian",NegTowers),0,1,NegTowers));//Median + //CHECK(METAlg::SubtractRho_MET(gTs,Form("RhoSubNeg%dMean",NegTowers),0,0,NegTowers));//Mean //Apply 3 sigma RMS - CHECK(METAlg::SubtractRho_MET(gTs,Form("RhoSubNeg%dMedianRMS",NegTowers),1,1,NegTowers)); - CHECK(METAlg::SubtractRho_MET(gTs,Form("RhoSubNeg%dMeanRMS",NegTowers),1,0,NegTowers)); + //CHECK(METAlg::SubtractRho_MET(gTs,Form("RhoSubNeg%dMedianRMS",NegTowers), 1, 1,NegTowers)); + //CHECK(METAlg::SubtractRho_MET(gTs,Form("RhoSubNeg%dMeanRMS",NegTowers), 1, 0,NegTowers)); CHECK(METAlg::Softkiller_MET(gTs, Form("SKNeg%d",NegTowers), NegTowers) ); //pileup subtracted SoftKiller (with avg rho) - CHECK(METAlg::JwoJ_MET(gTs,Form("JwoJNeg%d",NegTowers),m_pTcone_cut, NegTowers) ); //Jets without Jets - CHECK(METAlg::Pufit_MET(gTs,Form("PUfitNeg%d",NegTowers), NegTowers) ); //L1 version of PUfit, using gTowers + CHECK(METAlg::JwoJ_MET(gTs,Form("JwoJNeg%d",NegTowers),m_pTcone_cut, false, false, NegTowers) ); //Jets without Jets + //CHECK(METAlg::Pufit_MET(gTs,Form("PUfitNeg%d",NegTowers), NegTowers) ); //L1 version of PUfit, using gTowers } }//developer met else{ CHECK(METAlg::Baseline_MET(gTs,"gXENOISECUT",gT_noise, m_useNegTowers)); - CHECK(METAlg::JwoJ_MET(gTs,"gXEJWOJ",m_pTcone_cut,m_useNegTowers)); - CHECK(METAlg::SubtractRho_MET(gTs,"gXERHO",m_useRMS,m_useMedian,m_useNegTowers)); + CHECK(METAlg::JwoJ_MET(gTs,"gXEJWOJ",m_pTcone_cut,false, false, m_useNegTowers)); + //CHECK(METAlg::JwoJ_MET(gTs,"gXEJWOJRHO",m_pTcone_cut,false, true, m_useNegTowers)); CHECK(METAlg::Pufit_MET(gTs,"gXEPUFIT", m_useNegTowers) ); - CHECK(METAlg::MET_etaBins(gTs, "RhoSubEtaBins",true, true, false)); //simulating the 3 fpgas in the gFEX for rho subtraction + //CHECK(METAlg::SubtractRho_MET(gTs, "RhoSub_A", 1, true, false, false)); //simulating the 3 fpgas in the gFEX for rho subtraction + //CHECK(METAlg::SubtractRho_MET(gTs, "RhoSub_B", 2, true, false, false)); //simulating the 3 fpgas in the gFEX for rho subtraction + //CHECK(METAlg::SubtractRho_MET(gTs, "RhoSub_C", 3, true, false, false)); //simulating the 3 fpgas in the gFEX for rho subtraction + //CHECK(METAlg::Pufit_MET(gTs, "PUfit_A", 1, false)); //simulating the 3 fpgas in the gFEX for PUfit + //CHECK(METAlg::Pufit_MET(gTs, "PUfit_B", 2, false)); //simulating the 3 fpgas in the gFEX for PUfit + //CHECK(METAlg::Pufit_MET(gTs, "PUfit_C", 3,false)); //simulating the 3 fpgas in the gFEX for PUfit + CHECK(METAlg::SubtractRho_MET(gTs,"gXERHO", false, m_useRMS,false, m_useNegTowers)); + CHECK(MET_etaBins(gTs,"RhoSub_etaBins", true, true, false, false) ); + CHECK(MET_etaBins(gTs,"PUfit_etaBins", true, false, false, true) ); + CHECK(MET_etaBins(gTs,"JwoJ_etaBins", true, false, true, false) ); }//main definitions for simplicity return StatusCode::SUCCESS; } +StatusCode JGTowerReader::MET_etaBins(const xAOD::JGTowerContainer* gTs, TString metName, bool usegFEX, bool useRhoSub, bool useJwoJ, bool usePUfit){ + + //std::cout<<"Calculating MET in rings of eta..."< met = std::make_shared(); + + TString metA = ""; + TString metB = ""; + TString metC = ""; + + if(usegFEX){ + + + //xAOD::JGTowerContainer fpga_b; + //xAOD::JGTowerContainer fpga_c; + + ConstDataVector temp_a(SG::VIEW_ELEMENTS); + ConstDataVector temp_b(SG::VIEW_ELEMENTS); + ConstDataVector temp_c(SG::VIEW_ELEMENTS); + + //loop for splitting gTowers into 3 fpgas + for(unsigned int t = 0; t < gTs->size(); t++){ + const xAOD::JGTower* tower = gTs->at(t); + float eta = tower->eta(); + + //std::cout<<"tower eta = "<eta()< -2.5 && eta < 0){ + temp_a.push_back(tower); + } + if(eta > 0 && eta < 2.5){ + temp_b.push_back(tower); + } + if(TMath::Abs(eta) >= 2.5){ + temp_c.push_back(tower); + } + } + const xAOD::JGTowerContainer* fpga_a = temp_a.asDataVector(); + const xAOD::JGTowerContainer* fpga_b = temp_b.asDataVector(); + const xAOD::JGTowerContainer* fpga_c = temp_c.asDataVector(); + + if(fpga_a->size() == 0) std::cout<<"FPGA A is empty! "<size() == 0) std::cout<<"FPGA B is empty! "<size() == 0) std::cout<<"FPGA C is empty! "<mht_x; + float mst_ax = METAlg::m_METMap[metA]->mst_x; + float mht_ay = METAlg::m_METMap[metA]->mht_y; + float mst_ay = METAlg::m_METMap[metA]->mst_y; + + float mht_bx = METAlg::m_METMap[metB]->mht_x; + float mst_bx = METAlg::m_METMap[metB]->mst_x; + float mht_by = METAlg::m_METMap[metB]->mht_y; + float mst_by = METAlg::m_METMap[metB]->mst_y; + + float mht_cx = METAlg::m_METMap[metC]->mht_x; + float mst_cx = METAlg::m_METMap[metC]->mst_x; + float mht_cy = METAlg::m_METMap[metC]->mht_y; + float mst_cy = METAlg::m_METMap[metC]->mst_y; + + float mht_x = mht_ax + mht_bx + mht_cx; + float mht_y = mht_ay + mht_by + mht_cy; + float mst_x = mst_ax + mst_bx + mst_cx; + float mst_y = mst_ay + mst_by + mst_cy; + + //std::cout<<"MSTx = "<first.Data()))); - CHECK(evtStore()->record(METContAux,Form("%s_METAux.",it->first.Data()))); + CHECK(evtStore()->record(METContAux,Form("%s_METAux",it->first.Data()))); ATH_MSG_DEBUG("Recording EnergySumRoI with name " << it->first.Data() << "_MET"); } diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/src/METAlg.cxx b/Trigger/TrigT1/TrigT1CaloFexSim/src/METAlg.cxx index fdb5358ed6e..3384fdf479e 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/src/METAlg.cxx +++ b/Trigger/TrigT1/TrigT1CaloFexSim/src/METAlg.cxx @@ -15,64 +15,10 @@ #include "TrigT1CaloFexSim/Softkiller.h" #include "TrigT1CaloFexSim/JwoJ.h" #include "TrigT1CaloFexSim/Pufit.h" - std::map> METAlg::m_METMap; +std::map> METAlg::m_METMap; -//------------------------------------------------------------------------------------------------ -StatusCode METAlg::MET_etaBins(const xAOD::JGTowerContainer* towers, TString metName, bool usegFEX, bool useRhoSub, bool){ - float met_x = 0; - float met_y = 0; - std::shared_ptr met = std::make_shared(); - - if(usegFEX){ //if we are using the gFEX geometry - std::vector fpga_a; - std::vector fpga_b; - std::vector fpga_c; - - //loop for splitting gTowers into 3 fpgas - for(unsigned int t = 0; t < towers->size(); t++){ - const xAOD::JGTower* tower = towers->at(t); - float eta = tower->eta(); - if(eta > -2.4 && eta < 0){ - fpga_a.push_back(tower); - } - if(eta > 0 && eta < 2.4){ - fpga_b.push_back(tower); - } - if(TMath::Abs(eta) >= 2.4){ - fpga_c.push_back(tower); - } - } - - if(useRhoSub){ - std::vector met_a = rhoSub(fpga_a, false); - std::vector met_b = rhoSub(fpga_b, false); - std::vector met_c = rhoSub(fpga_c, false); - - float met_ax = met_a.at(0)*TMath::Cos(met_a.at(1)); - float met_ay = met_a.at(0)*TMath::Sin(met_a.at(1)); - - float met_bx = met_b.at(0)*TMath::Cos(met_b.at(1)); - float met_by = met_b.at(0)*TMath::Sin(met_b.at(1)); - - float met_cx = met_c.at(0)*TMath::Cos(met_c.at(1)); - float met_cy = met_c.at(0)*TMath::Sin(met_c.at(1)); - - met_x = met_ax + met_bx + met_cx; - met_y = met_ay + met_by + met_cy; - - float met_et = TMath::Sqrt(met_x*met_x + met_y*met_y); - float phi = TMath::ACos(met_x/met_et); - met->et = met_et; - met->phi = phi; - } - } - - if(m_METMap.find(metName)==m_METMap.end()) m_METMap[metName] = met; - return StatusCode::SUCCESS; -} //------------------------------------------------------------------------------------------------ StatusCode METAlg::Baseline_MET(const xAOD::JGTowerContainer*towers, TString metName, std::vector noise, bool useNegTowers){ - float met_x=0; float met_y=0; @@ -107,38 +53,83 @@ StatusCode METAlg::Baseline_MET(const xAOD::JGTowerContainer*towers, TString met return StatusCode::SUCCESS; } //------------------------------------------------------------------------------------------------ -StatusCode METAlg::SubtractRho_MET(const xAOD::JGTowerContainer* towers, TString metName, bool useRMS, bool useMedian, bool useNegTowers){ - +StatusCode METAlg::SubtractRho_MET(const xAOD::JGTowerContainer* towers, TString metName, bool useEtaBins, bool useRMS, bool useMedian, bool useNegTowers){ + float EtMiss = 0; + float n = 3.; float Ex = 0, Ey = 0, Ex_ = 0, Ey_ = 0; float threshold = 0; + float default_area = 0.04; + int fpga = 0; //can calculate rho as either the average or median gTower energy in the barrel - float rho = Rho_bar(towers,0); + if(metName.Contains("RhoSubA")) fpga = 1; + if(metName.Contains("RhoSubB")) fpga = 2; + if(metName.Contains("RhoSubC")) fpga = 3; + + //std::cout<size(); TH1F* h_Et = new TH1F("h_Et", "", 50, 0, 5000); if(useRMS){ + for(unsigned int t = 0; t < size; t++){ const xAOD::JGTower* tower = towers->at(t); + //float eta = tower->eta(); + float Et = tower->et(); + float eta = fabs(tower->eta()); + if(!useNegTowers && Et < 0) continue; h_Et->Fill(Et); } - threshold = h_Et->GetRMS(); + + threshold = n*(h_Et->GetRMS()); + //std::cout<<"3 sigma = "<at(t); + float Et = tower->et(); float phi = tower->phi(); - float eta = TMath::Abs(tower->eta()); + float eta = tower->eta(); + + //std::cout<<"Looking at tower ("<deta(); + //float dphi = tower->dphi(); + + float area = 0.;//deta*dphi/default_area; + if(eta >= -2.5 && eta < 0 && fpga == 1){ + if(eta > -2.4) area = 1.; + else area = 0.5; + } + if(eta > 0 && eta <= 2.5 && fpga == 2){ + if(eta < 2.4) area = 1.; + else area = 0.5; + } + if(fabs(eta) > 2.5 && fpga == 3){ + if(fabs(eta) < 3.2) area = 1.; + else area = 4.; + } + //std::cout<<"area = "< met = std::make_shared(); - met->phi=phi_met; + met->phi = phi_met; met->et = EtMiss; met->rho = rho; + + //std::cout<<"rho = "<size(); - - std::vector Et_values = Run_JwoJ(towers, pTcone_cut, useNegTowers); - + + float rho = 0; + if(useRho) rho = Rho_bar(towers, useEtaBins, 0, false); + else rho = 0; + std::vector Et_values = Run_JwoJ(towers, rho, pTcone_cut, useEtaBins, useNegTowers); + //set fit parameters for calculating MET //Look up table for parameters a,b,c depending on scalar sumEt //Optimized for resolution in external framework - float a; - float b; - float c; + float ax; + float bx; + float cx; + + float ay; + float by; + float cy; if(Et_values[0] <= 500*Gaudi::Units::GeV){ - a = 0.6; - b = 0.6; - c = 20.; + ax = 1.45; + bx = 1.15; + cx = 1.; + + ay = 1.45; + by = 1.1; + cy = 1.5; } else if(Et_values[0] <= 700*Gaudi::Units::GeV){ - a = 0.55; - b = 0.7; - c = 15.; + ax = 1.35; + bx = 1.05; + cx = -1.5; + + ay = 1.35; + by = 0.95; + cy = 0.; } else if(Et_values[0] <= 900*Gaudi::Units::GeV){ - a = 0.65; - b = 0.6; - c = 13.; + ax = 1.35; + bx = 1.0; + cx = 0.5; + + ay = 1.35; + by = 0.95; + cy = -1.; } else if(Et_values[0] <= 1100*Gaudi::Units::GeV){ - a = 0.75; - b = 0.55; - c = 10.; + ax = 1.3; + bx = 0.95; + cx = -1.; + + ay = 1.3; + by = 0.95; + cy = -1; } else if(Et_values[0] <= 1300*Gaudi::Units::GeV){ - a = 0.75; - b = 0.45; - c = 5.; + ax = 1.3; + bx = 0.9; + cx = 0.75; + + ay = 1.25; + by = 0.8; + cy = 0.5; } else if(Et_values[0] <= 1500*Gaudi::Units::GeV){ - a = 0.8; - b = 0.35; - c = 0.; + ax = 1.25; + bx = 0.8; + cx = 0.5; + + ay = 1.25; + by = 0.8; + cy = 0.5; } else if(Et_values[0] <= 1700*Gaudi::Units::GeV){ - a = 0.75; - b = 0.4; - c = -5.; + ax = 1.3; + bx = 0.75; + cx = 1.5; + + ay = 1.25; + by = 0.75; + cy = 0.5; } else{ - a = 0.7; - b = 0.45; - c = -10.; + ax = 1.25; + bx = 0.75; + cx = 0.75; + + ay = 1.25; + by = 0.75; + cy = 2.; } //a is hard term from gBlocks with pT > 25 GEV //b is total MET term computed from all gTowers - float EtMiss = a*Et_values[1] + b*Et_values[2] + c; + float Ex = ax*(Et_values[1])+ bx*Et_values[3] + cx; + float Ey = ay*(Et_values[2])+ by*Et_values[4] + cy; + + float EtMiss = TMath::Sqrt(Ex*Ex + Ey*Ey); + float phi = 0; + + float mht_phi = 0; + if(Et_values[5]>0) mht_phi = TMath::ACos(Et_values[1]/Et_values[5]); + if(Et_values[2] < 0) mht_phi = -mht_phi; + + float mst_phi = 0; + if(Et_values[6]>0) mst_phi = TMath::ACos(Et_values[3]/Et_values[6]); + if(Et_values[4] < 0) mst_phi = -mst_phi; + + if(EtMiss > 0) phi = TMath::ACos(Ex/EtMiss); + + if(Ey < 0) phi = -phi; std::shared_ptr met = std::make_shared(); - met->phi=0; + met->phi=phi; met->et =EtMiss; - met->rho = 0; - met->mst = Et_values[2]; - met->mht = Et_values[1]; + //met->rho = rho; + //std::cout<<"Hard term = "<eta(); float Et = tower->et(); - float eta = fabs(tower->eta()); if(!useNegTowers && Et < 0) continue; h_Et->Fill(Et); @@ -405,7 +403,6 @@ float METAlg::Rho_avg_etaRings(const xAOD::JGTowerContainer* towers, int fpga, b const xAOD::JGTower* tower = towers->at(i); float eta = tower->eta(); - float phi = tower->phi(); float Et = tower->et(); if(!useNegTowers && Et < 0) continue; -- GitLab From 25e697c6040699c1c9c6e36390b11621422e2e81 Mon Sep 17 00:00:00 2001 From: amyers Date: Mon, 24 Jun 2019 13:32:38 +0200 Subject: [PATCH 3/5] Saved gCaloTowers and pileup subtracted towers to SG, fixed logic issues --- Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/Rho.h | 6 +++--- Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx | 8 +++++++- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/Rho.h b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/Rho.h index 5e46552ce77..ee30a461073 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/Rho.h +++ b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/Rho.h @@ -58,7 +58,7 @@ float Rho_bar(const xAOD::JGTowerContainer* towers, const bool useEtaBins, int f for(unsigned int i = 0; i < size; i++){ const xAOD::JGTower* tower = towers->at(i); - float eta = TMath::Abs(tower->eta()); + float eta = tower->eta(); float Et = tower->et(); if(!useNegTowers && Et < 0) continue; @@ -72,7 +72,7 @@ float Rho_bar(const xAOD::JGTowerContainer* towers, const bool useEtaBins, int f else area_b += 0.5; } if(fabs(eta) > 2.5){ - if(eta < 3.2) area_c += 1.; + if(fabs(eta) < 3.2) area_c += 1.; else area_c += 4.; } @@ -82,7 +82,7 @@ float Rho_bar(const xAOD::JGTowerContainer* towers, const bool useEtaBins, int f //length++; } }else{ - if(eta < 2.4){ + if(fabs(eta) < 2.4){ if(Et < et_max){ rho += Et; length++; diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx b/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx index da2bbd8cac0..fd53555ac25 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx +++ b/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx @@ -211,6 +211,7 @@ StatusCode JGTowerReader::execute() { xAOD::JGTowerAuxContainer* gCaloTowersAux = new xAOD::JGTowerAuxContainer(); xAOD::JGTowerContainer* gCaloTowers = new xAOD::JGTowerContainer(); gCaloTowers->setStore(gCaloTowersAux); + for(unsigned t = 0; t < gTowers->size(); t++){ const xAOD::JGTower* gt_em = gTowers->at(t); const float eta = gt_em->eta(); @@ -331,6 +332,10 @@ StatusCode JGTowerReader::execute() { ATH_MSG_DEBUG ("Algorithm passed"); + //manage containers that have been created: save gCaloTowers and pu_sub to SG + CHECK(evtStore()->record(gCaloTowers,"gCaloTowers")); + CHECK(evtStore()->record(pu_sub,"pu_subTowers")); + setFilterPassed(true); //if got here, assume that means algorithm passed return StatusCode::SUCCESS; } @@ -703,6 +708,7 @@ StatusCode JGTowerReader::MET_etaBins(const xAOD::JGTowerContainer* gTs, TString } if(METAlg::m_METMap.find(metName) == METAlg::m_METMap.end()) METAlg::m_METMap[metName] = met; + return StatusCode::SUCCESS; } @@ -763,7 +769,7 @@ StatusCode JGTowerReader::ProcessObjects(){ (*acc_mht)(*METCont) = met->mht; (*acc_mst)(*METCont) = met->mst; CHECK(evtStore()->record(METCont,Form("%s_MET",it->first.Data()))); - CHECK(evtStore()->record(METContAux,Form("%s_METAux",it->first.Data()))); + CHECK(evtStore()->record(METContAux,Form("%s_METAux.",it->first.Data()))); ATH_MSG_DEBUG("Recording EnergySumRoI with name " << it->first.Data() << "_MET"); } -- GitLab From 9d452566eee50d782e90ecd1a934506ab428a167 Mon Sep 17 00:00:00 2001 From: amyers Date: Mon, 24 Jun 2019 13:45:13 +0200 Subject: [PATCH 4/5] Saved Aux's to SG --- Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx b/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx index fd53555ac25..9386ede4126 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx +++ b/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx @@ -334,7 +334,9 @@ StatusCode JGTowerReader::execute() { //manage containers that have been created: save gCaloTowers and pu_sub to SG CHECK(evtStore()->record(gCaloTowers,"gCaloTowers")); + CHECK(evtStore()->record(gCaloTowersAux, "gCaloTowersAux.")); CHECK(evtStore()->record(pu_sub,"pu_subTowers")); + CHECK(evtStore()->record(pu_subAux, "pu_subTowersAux.")); setFilterPassed(true); //if got here, assume that means algorithm passed return StatusCode::SUCCESS; -- GitLab From a2d2100f1c76e4578197f21461ea61091375d16c Mon Sep 17 00:00:00 2001 From: amyers Date: Mon, 24 Jun 2019 15:21:27 +0200 Subject: [PATCH 5/5] fixed fabs(eta) statement --- Trigger/TrigT1/TrigT1CaloFexSim/src/METAlg.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/src/METAlg.cxx b/Trigger/TrigT1/TrigT1CaloFexSim/src/METAlg.cxx index bad3a236d8d..48b4b3c5cfe 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/src/METAlg.cxx +++ b/Trigger/TrigT1/TrigT1CaloFexSim/src/METAlg.cxx @@ -123,7 +123,7 @@ StatusCode METAlg::SubtractRho_MET(const xAOD::JGTowerContainer* towers, TString if(useEtaBins) Et_sub = TMath::Abs(Et) - area*rho; else{ - if(eta < 2.4) Et_sub = TMath::Abs(Et) - rho; + if(fabs(eta < 2.4)) Et_sub = TMath::Abs(Et) - rho; else Et_sub = TMath::Abs(Et) - 4*rho; } -- GitLab