Merge branch 'MET_algs' into '21.3'

Implemented gFEX MET algorithms: Rho subtraction, Softkiller, Jets without Jets,…

See merge request atlas/athena!15835

Former-commit-id: ba887415edc08b9fab2f598bc50c1e2fef129a7e

Trigger/TrigT1/TrigT1CaloFexSim/CMakeLists.txt

@@ -28,7 +28,8 @@ atlas_depends_on_subdirs( PUBLIC
                           Trigger/TrigAnalysis/TrigAnalysisInterfaces)
 
 # External dependencies:
-find_package( ROOT COMPONENTS Core Tree MathCore Hist RIO pthread )
+find_package( ROOT COMPONENTS Core Tree MathCore Hist Matrix RIO pthread)
+
 
 # Component(s) in the package:
 atlas_add_component( TrigT1CaloFexSim

Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JGTowerReader.h

@@ -42,6 +42,7 @@ class JGTowerReader: public ::AthAlgorithm {
   virtual StatusCode  HistBookFill(const TString name, Int_t nbinsx, Double_t xbin_down, Double_t xbin_up, float xvalue, float wei);
   virtual StatusCode  HistBookFill(const TString name, Int_t nbinsx, const Double_t* xbins, float xvalue,float wei);
  private: 
+  bool m_vetoBCID;
   bool m_outputNoise;
   bool m_debugJetAlg;
   bool m_dumpTowersEtaPhi;
@@ -66,6 +67,15 @@ class JGTowerReader: public ::AthAlgorithm {
   float m_gJet_r;
   std::string m_noise_file;
  
+  //job options for gFEX MET algorithms
+  bool m_useRMS;
+  bool m_useMedian;
+  bool m_useNegTowers;
+  bool m_combine_rhoNoise;
+  bool m_combine_skNoise;
+  bool m_combine_jwojNoise;
+  float m_pTcone_cut;
+
   const CaloCell_SuperCell_ID* m_scid;
   const JTower_ID* m_jTowerId;
   const GTower_ID* m_gTowerId;
@@ -85,7 +95,22 @@ class JGTowerReader: public ::AthAlgorithm {
   JetAlg::Seed*   jJetSeeds = new JetAlg::Seed;
   JetAlg::Seed*   gSeeds=new JetAlg::Seed;
   METAlg::MET*    jMET=new  METAlg::MET;
-  METAlg::MET*    gMET=new  METAlg::MET;
+
+  std::map<TString, METAlg::MET*> met_algs;
+
+  METAlg::MET*    gMET = new  METAlg::MET;
+  METAlg::MET*    gMET_rho = new  METAlg::MET;
+  METAlg::MET*    gMET_sk = new  METAlg::MET;
+  METAlg::MET*    gMET_jwoj = new  METAlg::MET;
+  METAlg::MET*    gMET_pufit = new  METAlg::MET;
+
+  /* met_algs["rho_sub"] = gMET_rho;  //pileup subtraction
+  met_algs["SK"] = gMET_sk;      //softkiller
+  met_algs["JwoJ"] = gMET_jwoj;    //jets without jets
+  met_algs["Noise"] = gMET;  //4 sigma noise cut
+  met_algs["PUfit"] = gMET_pufit;  //pufit
+*/  
+
   std::vector<JetAlg::L1Jet>  jL1Jets;
   std::vector<JetAlg::L1Jet>  jJet_L1Jets;
   std::vector<JetAlg::L1Jet>  gL1Jets;

Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JetAlg.h

@@ -60,7 +60,8 @@ class JetAlg{
   static StatusCode SeedGrid(const xAOD::JGTowerContainer*towers, JetAlg::Seed*seeds, bool &m_dumpSeedsEtaPhi);
   static StatusCode SeedFinding(const xAOD::JGTowerContainer*towers, JetAlg::Seed*seeds, float seed_size,float range, std::vector<float> noise, bool m_debug = false);
 
-  static StatusCode BuildJet(const xAOD::JGTowerContainer*towers, JetAlg::Seed*seeds, std::vector<JetAlg::L1Jet> &js, float jet_size, std::vector<float> noise, bool m_debug = false);
+  static StatusCode BuildFatJet(const xAOD::JGTowerContainer towers, std::vector<JetAlg::L1Jet>& js, float jet_r, std::vector<float> noise); 
+  static StatusCode BuildJet(const xAOD::JGTowerContainer*towers, JetAlg::Seed*seeds, std::vector<JetAlg::L1Jet> &js, float jet_size, std::vector<float> noise, bool m_debug);
   static StatusCode BuildRoundJet(const xAOD::JGTowerContainer*towers, JetAlg::Seed*seeds, std::vector<JetAlg::L1Jet> &js, float jet_size, std::vector<float> noise, bool m_debug = false);
 
 };

Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JwoJ.h

+/**
+
+ *   Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+
+ */
+
+#ifndef TRIGT1CALOFEXSIM_JWOJ_H
+#define TRIGT1CALOFEXSIM_JWOJ_H
+/*
+ * Class   : JwoJ
+ *
+ * Author  : Myers, Ava (amyers@cern.ch)
+ *
+ * Date    : Oct 2018 - Created class
+ *
+ * Implements the Jets without Jets algorithm for MET reconstruction. Optimized for resolution (Rebecca Linck): Uses total reconstructed MET as the "soft term". Includes all helper methods for MET calculation using JwoJ
+ */
+#include <algorithm>
+#include "TMath.h"
+#include "Objects.h"
+#include "GaudiKernel/SystemOfUnits.h"
+
+
+/**
+ *@brief Takes JGTowers from a JGTowerContainer and combines them into gBlocks. Used for threshold on the "cone" (gBlock) pt
+ *@return @c void 
+ */
+void BuildBlocksFromTowers(std::vector<TowerObject::Block>& blocks, const xAOD::JGTowerContainer towers, const int blockRows, const int blockCols, const bool useNegTowers){
+  
+  blocks.clear();
+  TowerObject::TowerGrid grid = TowerObject::TowerGrid(towers);
+  
+  for(const xAOD::JGTower* seed: towers){
+    int seedIndex = std::find(towers.begin(), towers.end(), seed) - towers.begin();
+    
+    std::vector<int> 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);
+    double neighbor_pt = 0;
+    
+    for(const int& neighborIndex: neighbors){
+      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;
+    }
+    
+    TowerObject::Block block(block_pt, seed->eta(), seed->phi(), 0.0);
+    block.seedIndex(seedIndex);
+    block.numEta(blockCols);
+    block.numPhi(blockRows);
+    block.area(block_area);
+    block.numConstituents(neighbors.size());
+    
+    blocks.push_back(block);
+  }
+  
+  std::sort(blocks.rbegin(), blocks.rend());
+}
+
+/**
+ *@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<float>
+ */
+std::vector<float> Run_JwoJ(const xAOD::JGTowerContainer* towers, float pTcone_cut, bool useNegTowers){
+  
+  std::vector<TowerObject::Block> blocks;
+  BuildBlocksFromTowers(blocks, *towers, 3, 3, useNegTowers);
+  
+  pTcone_cut*=Gaudi::Units::GeV;
+  std::vector<float> Et_vals;
+  float Ht = 0; float Htx = 0; float Hty = 0;
+  float Et = 0; float Etx = 0; float Ety = 0;
+  float Et_tot = 0; float Etx_tot = 0; float Ety_tot = 0;
+  
+  float scalar_Et = 0;
+  
+  for(unsigned int b = 0; b < blocks.size(); b++){
+    float block_phi = blocks[b].Phi();
+    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);
+    
+    if(TMath::Abs(blocks[b].Eta()) < 2.4){
+      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{
+      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;
+    Etx_tot += block_etx;
+    Ety_tot += block_ety;
+  }
+  float MHT = TMath::Sqrt(Htx*Htx + Hty*Hty);
+  float MET = TMath::Sqrt(Etx*Etx + Ety*Ety);
+  Et_tot = TMath::Sqrt(Etx_tot*Etx_tot + Ety_tot*Ety_tot);
+
+  //convert to GeV for fitting process  
+  Et_vals.push_back(scalar_Et);
+  Et_vals.push_back(MHT);
+  Et_vals.push_back(MET);
+  Et_vals.push_back(Et_tot);
+  
+  return Et_vals;
+}
+
+#endif
+

Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/METAlg.h

@@ -2,6 +2,22 @@
  *   Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
 */
 
+#ifndef TRIGT1CALOFEXSIM_METALG_H
+#define TRIGT1CALOFEXSIM_METALG_H
+
+/*
+ * Class    : METAlg
+ *
+ * Authors  : Lin, Chiao-Ying (cylin@cern.ch); Myers, Ava (amyers@cern.ch)
+ *
+ * Date     : Oct 2018 - Updated class
+ *
+ * Includes algorithms for MET reconstruction and pileup suppression.
+ */
+
+#include "GaudiKernel/MsgStream.h"
+#include "AthenaBaseComps/AthAlgTool.h"
+
 #include "CaloIdentifier/GTower_ID.h"
 #include "CaloEvent/CaloCellContainer.h"
 #include "xAODTrigL1Calo/JGTower.h"
@@ -15,16 +31,35 @@
 #include "TH1.h"
 #include "TH2.h"
 
-class METAlg{
 
+class METAlg{
+  
  public:
   struct MET{
     float phi;
     float et;
-
+    
   };
-
-
-  static StatusCode BuildMET(const xAOD::JGTowerContainer*towers,METAlg::MET* met,std::vector<float> noise);
-
+  /**
+   *@brief Calculate MET using a fixed 4 sigma noise cut
+   */
+  static StatusCode BuildMET(const xAOD::JGTowerContainer*towers,METAlg::MET* met,std::vector<float> noise, bool useNegTowers);
+  /**
+   *@brief Calculates MET with pileup subtraction
+   */
+  static StatusCode SubtractRho_MET(const xAOD::JGTowerContainer* towers, METAlg::MET* met, bool useRMS, bool useMedian, bool useNegTowers);
+  /**
+   *@brief Calculates MET with Softkiller
+   */
+  static StatusCode Softkiller_MET(const xAOD::JGTowerContainer* towers, METAlg::MET* met, bool useNegTowers);
+  /**
+   *@brief Calculates MET with Jets without Jets
+   */
+  static StatusCode JwoJ_MET(const xAOD::JGTowerContainer* towers, METAlg::MET* met, float pTcone_cut, bool useNegTowers);
+  /**
+   *@brief Calculates MET using PUfit
+   */
+  static StatusCode Pufit_MET(const xAOD::JGTowerContainer* towers, METAlg::MET* met, bool useNegTowers);
 };
+
+#endif

Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/Objects.h

+/*
+
+ *   Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+
+Objects.h: Contains definitions of objects necessary for MET algorithms
+
+Author: Ava Myers (amyers@cern.ch)
+
+Objects: Object: template four momentum object holding pt, eta, phi, m for definitions of more complex objects
+
+         Block: g(j)Block required for Jets without Jets algorithm. Designed to be 3x3 block around a seed tower
+
+         TowerGrid: Used by Block object. Defines coordinates/indices of towers and neighbors to the seed tower
+
+	 Patch: Tower "patch" required for PUfit algorithm. A simpler implementation of Block/TowerGrid for quick use. Called a tower in PUfit algorithm, but renaming as Patch to avoid confusion
+ */
+
+#include "Math/Vector4D.h"
+#include "TVector2.h"
+#include "math.h"
+#include "TMath.h"
+#include <vector>
+#include <algorithm>
+#include "TrigT1CaloFexSim/JGTower.h"
+
+namespace TowerObject{
+  typedef ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<float>> FourMom_t;
+  
+  //primarily to make the interface consistent
+  class Object: public FourMom_t{
+  private:
+    float m_area = 0.0;
+    
+  public:
+  Object(): FourMom_t(0.0, 0.0, 0.0, 0.0){};
+  Object(float pt, float eta, float phi, float m):
+    FourMom_t(pt, eta, phi, m){};
+    ~Object(){};
+    
+    float pt() const {return this->Pt();};
+    void pt(float pt) {this->SetPt(pt);};
+    float eta() const {return this->Eta();};
+    void eta(float eta) {this->SetEta(eta);};
+    float phi() const {return this->Phi();};
+    void phi(float phi) {this->SetPhi(phi);};
+    float m() const {return this->M();};
+    void m(float m) {this->SetM(m);};
+    float area() const {return this->M();};
+    void area(float a) {m_area = a;};
+    
+    // comparison operators for easy sorting                                                                                    
+    friend bool operator <  (const Object &o1, const Object &o2){ return o1.pt() <  o2.pt(); };
+    friend bool operator >  (const Object &o1, const Object &o2){ return o1.pt() >  o2.pt(); };
+    friend bool operator<= (const Object &o1, const Object &o2){ return o1.pt() <= o2.pt(); };
+    friend bool operator>= (const Object &o1, const Object &o2){ return o1.pt() >= o2.pt(); };
+    
+    float deltaR(const Object& other) const{
+      float deta = this->deltaEta(other);
+      float dphi = this->deltaPhi(other);
+      return sqrt(deta*deta + dphi*dphi);
+    };
+    
+    float deltaEta(const Object& other) const{
+      return this->eta() - other.eta();
+    };
+    float deltaPhi(const Object& other) const{
+      return TVector2::Phi_mpi_pi(this->phi() - other.phi());
+    };
+  
+  };
+
+  class Block : public Object{
+  private: 
+    int m_seedIndex = -1;
+    int m_numEta;
+    int m_numPhi;
+    int m_numConstituents = 0;
+    
+  public:
+    //Block: Object(){};
+  Block(float pt, float eta, float phi, float m): Object(pt, eta, phi, m){};
+    
+    float pt_cal() const {return 2700 + (1./0.8189)*this->Object::pt();}; //values have been previously derived for gFEX calibrations
+    //    float pt(bool calib = false) const {return (calibration)?pt_cal():this->Object::pt();}
+    int seedIndex() const {return m_seedIndex;};
+    void seedIndex(int index) {m_seedIndex = index;};
+    int numEta() const {return m_numEta;};
+    void numEta(int n) {m_numEta = n;};
+    int numPhi() const {return m_numPhi;};
+    void numPhi(int n) {m_numPhi = n;};
+     int numConstituents() const {return m_numConstituents;};
+    void numConstituents(int n) {m_numConstituents = n;};
+  };
+  
+  class TowerGrid{
+  public:
+    //build the grid
+    TowerGrid(const xAOD::JGTowerContainer& towers){
+      grid = std::vector<int>(grid_eta*grid_phi, -1);
+      for(const xAOD::JGTower* tower: towers){
+	int towerIndex = std::find(towers.begin(), towers.end(), tower) - towers.begin();
+	this->add_tower(*tower, towerIndex);
+      }
+    };
+    
+    void add_tower(const xAOD::JGTower& tower, int towerIndex){
+      int gridIndex = get_index(tower);
+      if(gridIndex >= 0){
+	grid[gridIndex] = towerIndex;
+      }
+    };
+    
+    //normalization functions
+    std::pair<int, int> normalize_coord(std::pair<int, int> coord) const{
+      while(coord.second < 0) coord.second += grid_phi;
+      return coord;
+    };
+    
+    int normalize_index(int index) const{
+      while(index < 0) index += grid_eta*grid_phi;
+      return index;
+    };
+  
+    //get index functions
+    int get_index(const xAOD::JGTower& tower) const{
+      return get_index(get_coord(tower));
+    };
+    int get_index(std::pair<int, int> coord) const{
+      coord = normalize_coord(coord);
+      if(coord.first < 0 || coord.first >= grid_eta) return -1;
+      return coord.second*grid_eta + coord.first;
+    };
+    
+    //get coordinate functions
+    std::pair<int, int> get_coord(const xAOD::JGTower& tower) const{
+      int ieta(-1);
+      for(const float& eta: eta_bins){
+	if(tower.eta() < eta) break;
+	ieta++;
+      }
+      int iphi(-1);
+      for(const float& phi : phi_bins){
+	if(tower.phi() < phi) break;
+	iphi++;
+      }
+      return std::pair<int, int>(ieta, iphi);
+    };
+    
+    std::pair<int, int> get_coord(int index) const{
+      index = normalize_index(index);
+      int iphi = std::floor(index/grid_eta);
+      int ieta = index - iphi*grid_eta;
+      return std::pair<int, int>(ieta, iphi);
+    };
+    
+    //neighbors functions
+    std::vector<int> neighbors(const xAOD::JGTower& tower, int nrows = 3, int ncols = 3) const{
+      const std::pair<int, int> coord = get_coord(tower);
+      return neighbors(coord, nrows, ncols);
+    };
+    
+    std::vector<int> neighbors(const std::pair<int, int> coord, int nrows, int ncols) const{
+      std::vector<int> neighbors;
+      if(nrows%2 == 0 || ncols%2 == 0) return neighbors;
+      
+      for(int irow = -(nrows-1)/2; irow <= (nrows-1)/2; irow++){
+	for(int icol = -(ncols-1)/2; icol <= (ncols-1)/2; icol++){
+	  //don't include the target tower itself
+	  if(irow == 0 && icol == 0) continue;
+	  std::pair<int, int> neighbor(coord.first + icol, coord.second + irow);
+	  int index = this->at(normalize_coord(neighbor));
+	  //only if the tower exists in the grid
+	  if(index >= 0){
+	    neighbors.push_back(index);
+	  }
+	}
+      }
+      return neighbors;
+    };
+    
+    //at functions
+    int at(const std::pair<int, int>& coord) const{
+      return this->at(get_index(coord));
+    };
+    int at(int index) const{
+      if(index < 0) return -1;
+      if(index >= static_cast<int>(grid.size())) return -1;
+      return grid.at(index);
+    };
+    
+  private:
+    const int grid_eta = 24;
+    const int grid_phi = 32;
+    std::vector<int> grid;
+    const std::vector<float> eta_bins = {
+      -2.4000,
+      -2.2000,
+      -2.0000,
+      -1.8000,
+      -1.6000,
+      -1.4000,
+      -1.2000,
+      -1.0000,
+      -0.8000,
+      -0.6000,
+      -0.4000,
+      -0.2000,
+      0.0000,
+      0.2000,
+      0.4000,
+      0.6000,
+      0.8000,
+      1.0000,
+      1.2000,
+      1.4000,
+      1.6000,
+      1.8000,
+      2.0000,
+      2.2000,
+      2.4000
+    };
+    
+    const std::vector<float> phi_bins = {
+      -3.1416,
+      -2.9452,
+      -2.7489,
+      -2.5525,
+      -2.3562,
+      -2.1598,
+      -1.9635,
+      -1.7671,
+      -1.5708,
+      -1.3744,
+      -1.1781,
+      -0.9817,
+      -0.7854,
+      -0.5890,
+      -0.3927,
+      -0.1963,
+      0.0000,
+      0.1963,
+      0.3927,
+      0.5890,
+      0.7854,
+      0.9817,
+      1.1781,
+      1.3744,
+      1.5708,
+      1.7671,
+      1.9635,
+      2.1598,
+      2.3562,
+      2.5525,
+      2.7489,
+      2.9452,
+      3.1416
+    };
+  }; 
+}

Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/Pufit.h

+/*
+ *   Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef TRIGT1CALOFEXSIM_PUFIT_H
+#define TRIGT1CALOFEXSIM_PUFIT_H
+/*
+ * Class   : Pufit
+ *
+ * Author  : Myers, Ava (amyers@cern.ch)
+ *
+ *Implements the Pufit algorithm for MET reconstruction and pileup subtraction. Includes a helper class for gTower patches
+ */
+#include "CaloIdentifier/GTower_ID.h"
+#include "CaloEvent/CaloCellContainer.h"
+#include "xAODTrigL1Calo/JGTower.h"
+#include "xAODTrigL1Calo/JGTowerContainer.h"
+#include "xAODTrigL1Calo/JGTowerAuxContainer.h"
+#include "xAODEventInfo/EventInfo.h"
+#include "JetCalibTools/IJetCalibrationTool.h"
+#include "JetInterface/IJetUpdateJvt.h"
+#include "TrigAnalysisInterfaces/IBunchCrossingTool.h"
+#include "Identifier/IdentifierHash.h"
+
+#include <TMatrixD.h>
+#include <TMath.h>
+#include <TLorentzVector.h>
+#include <vector>
+
+/**
+ *@brief Helper class for patches of towers
+ */
+struct Patch{
+  float sumEt{0.};
+  float px{0.};
+  float py{0.};
+  float pt() {return sqrt(px*px + py*py);}
+  float cosPhi() {return pt() > 0 ? px/pt() : 1.;}
+  float sinPhi() {return pt() > 0 ? py/pt() : 0.;}
+  Patch & add(float Et, float Phi){
+    float sinPhi = TMath::Sin(Phi);
+    float cosPhi = TMath::Cos(Phi);
+    sumEt += Et;
+    px += Et*cosPhi;
+    py += Et*sinPhi;
+    return *this;
+  };
+};
+
+/**
+ *@brief Namespace for "global" variables in PUfit
+ */
+namespace PUfitVar{
+  float m_targetTowerWidth = 0.7;
+  float m_maxEta = 5.;
+  float m_caloResSqrtTerm = 15.81;
+  float m_caloResFloor = 50;
+  float m_nSigma = 3.0;