diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JGTowerReader.h b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JGTowerReader.h index db4343368ca84ccb70317d805c34fae6eeda40a7..e2ccd7d12a78c186d01c8145a06639b999a41de4 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 bb96f74ccef5e4f5f256d6854b0a62e7c5d6cee6..162e23b52c0c301d846a6b24342d138508fb2d6c 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JwoJ.h +++ b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JwoJ.h @@ -35,6 +35,7 @@ void BuildBlocksFromTowers(std::vector& blocks, const xAOD:: std::vector neighbors = grid.neighbors(*seed, blockRows, blockCols); float seed_Et = seed->et(); + if(!useNegTowers) seed_Et = TMath::Abs(seed_Et); double block_area(0.0); double block_pt(seed_Et); @@ -44,6 +45,7 @@ void BuildBlocksFromTowers(std::vector& blocks, const xAOD:: const xAOD::JGTower* neighbor = towers.at(neighborIndex); block_area += neighbor->deta()*neighbor->dphi(); neighbor_pt = neighbor->et(); + if(!useNegTowers) neighbor_pt = TMath::Abs(neighbor_pt); block_pt += neighbor_pt; } @@ -65,7 +67,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 +85,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 +101,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 +144,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 f1f2ecca895800ecfb4bef409ab03bbe2b38310e..15bce7c7ff2065ef57621080a55b5722ae05e458 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 a8a07d9068709036ceae3f1e44c95aae35dc4f1a..ee30a4610738dca39e7453d5f5730080ef3f6061 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 eta = 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(fabs(eta) < 3.2) area_c += 1.; + else area_c += 4.; + } + + if(useEtaBins){ if(Et < et_max){ rho += Et; - length++; + //length++; + } + }else{ + if(fabs(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 6655fd6740fb8b37f9bdf571b4865df2cbf8c837..9386ede4126c43e64c88cf485956c0bb8026686e 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,108 @@ 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 +323,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 @@ -219,6 +332,12 @@ 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(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; } @@ -335,29 +454,265 @@ 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 = "<> 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,81 @@ 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; + 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(); + 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 = "<