ATLAS Offline Software
AnomDetVAE.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 /*********************************
5  * AnomDetVAE.cxx
6  * Created by Sagar Addepalli on 25/11/2024.
7  *
8  * @brief algorithm uses a variational auto-encoder based anomaly detection network
9  * Uses 6 jets, 4 muons, 4 taus, and MET to calculate the anomaly score
10  * based on the KL-divergence in a lower dimension latent space
11  *
12  * @param ScaleSqr1
13  * @param ScaleSqr2
14  * @param ScaleSqr3
15  * @param AnomalyScoreThresh
16 **********************************/
17 
18 #include <cmath>
19 
20 #include "L1TopoAlgorithms/AnomDetVAE.h"
21 #include "L1TopoCommon/Exception.h"
22 #include "L1TopoInterfaces/Decision.h"
23 #include "L1TopoSimulationUtils/Helpers.h"
24 #include <VAENetwork.h>
25 
26 REGISTER_ALG_TCS(ADVAE_2A)
27 
28 
29 TCS::ADVAE_2A::ADVAE_2A(const std::string & name) : DecisionAlg(name)
30 {
31  defineParameter("InputWidth1", 6);
32  defineParameter("InputWidth2", 6);
33  defineParameter("InputWidth3", 6);
34  defineParameter("InputWidth4", 1);
35  defineParameter("MaxTob1", 6);
36  defineParameter("MaxTob2", 4);
37  defineParameter("MaxTob3", 4);
38  defineParameter("MaxTob4", 1);
39  defineParameter("NumResultBits", 2);
40 
41  // Version parameter, used for L1TopoFW book-keeping, no practical application
42  defineParameter("ADVAEVersion", 1);
43  // minEt cuts, one per input list (TOBs failing these are set to ET = eta = phi = 0)
44  defineParameter("MinET1",0);
45  defineParameter("MinET2",0);
46  defineParameter("MinET3",0);
47  defineParameter("MinET4",0);
48  // Configurable scale for the mu parameters in the latent space and AD score threshold
49  // The value used is sum of squares of the NN result vector elements, bit-shifted by 8 to get all decimal bits
50  defineParameter("ScaleSqr1",128,0);
51  defineParameter("ScaleSqr2",128,0);
52  defineParameter("ScaleSqr3",128,0);
53  defineParameter("AnomalyScoreThresh", 1000000, 0);
54  defineParameter("ScaleSqr1",128,1);
55  defineParameter("ScaleSqr2",128,1);
56  defineParameter("ScaleSqr3",128,1);
57  defineParameter("AnomalyScoreThresh", 1000000, 1);
58 
59  setNumberOutputBits(2);
60 }
61 
62 TCS::ADVAE_2A::~ADVAE_2A(){}
63 
64 
65 TCS::StatusCode
66 TCS::ADVAE_2A::initialize() {
67  p_NumberLeading1 = parameter("InputWidth1").value();
68  p_NumberLeading2 = parameter("InputWidth2").value();
69  p_NumberLeading3 = parameter("InputWidth3").value();
70  p_NumberLeading4 = parameter("InputWidth4").value();
71 
72  if(parameter("MaxTob1").value() > 0) p_NumberLeading1 = parameter("MaxTob1").value();
73  if(parameter("MaxTob2").value() > 0) p_NumberLeading2 = parameter("MaxTob2").value();
74  if(parameter("MaxTob3").value() > 0) p_NumberLeading3 = parameter("MaxTob3").value();
75  if(parameter("MaxTob4").value() > 0) p_NumberLeading4 = parameter("MaxTob4").value();
76 
77  p_minEt1 = parameter("MinET1").value();
78  p_minEt2 = parameter("MinET2").value();
79  p_minEt3 = parameter("MinET3").value();
80  p_minEt4 = parameter("MinET4").value();
81 
82  for(unsigned int i=0; i<numberOutputBits(); ++i) {
83  p_ScaleSqr1[i] = parameter("ScaleSqr1", i).value();
84  p_ScaleSqr2[i] = parameter("ScaleSqr2", i).value();
85  p_ScaleSqr3[i] = parameter("ScaleSqr3", i).value();
86  p_AnomalyScoreThresh[i] = parameter("AnomalyScoreThresh", i).value();
87  }
88 
89  TRG_MSG_INFO("number output : " << numberOutputBits());
90 
91  // book histograms
92  for(unsigned int i=0; i<numberOutputBits(); ++i) {
93  std::string hname_accept = "hAnomalyScore_accept_bit"+std::to_string((int)i);
94  std::string hname_reject = "hAnomalyScore_reject_bit"+std::to_string((int)i);
95  // score
96  bookHist(m_histAccept, hname_accept, "ADScore", 2000, 0, 2000000);
97  bookHist(m_histReject, hname_reject, "ADScore", 2000, 0, 2000000);
98  }
99 
100  return StatusCode::SUCCESS;
101 }
102 
103 
104 TCS::StatusCode
105 TCS::ADVAE_2A::processBitCorrect( const std::vector<TCS::TOBArray const *> & input,
106  const std::vector<TCS::TOBArray *> & output,
107  Decision & decision )
108 {
109 
110 
111  if( input.size() == 4) {
112 
113  TCS::TOBArray const* jets = input[0];
114  TCS::TOBArray const* taus = input[1];
115  TCS::TOBArray const* mus = input[2];
116  TCS::TOBArray const* met = input[3];
117  TRG_MSG_DEBUG("Number of jets are " << (*jets).size());
118  TRG_MSG_DEBUG("Number of taus are " << (*taus).size());
119  TRG_MSG_DEBUG("Number of mus are " << (*mus).size());
120  TRG_MSG_DEBUG("Number of met are " << (*met).size());
121 
122  //check for ambiguous sorting and set corresponding flag if an ambiguity is found
123  bool hasAmbiguousInputs = TSU::isAmbiguousAnywhere(jets, p_NumberLeading1, p_minEt1)
124  || TSU::isAmbiguousAnywhere(taus, p_NumberLeading2, p_minEt2)
125  || TSU::isAmbiguousAnywhere(mus, p_NumberLeading3, p_minEt3)
126  || TSU::isAmbiguousAnywhere(met, p_NumberLeading4, p_minEt4);
127 
128  std::vector<u_int> jet_pt(6,0), tau_pt(4,0), mu_pt(4,0), met_pt(1,0);
129  std::vector<int> jet_eta(6,0), tau_eta(4,0), mu_eta(4,0); //no met_eta
130  std::vector<int> jet_phi(6,0), tau_phi(4,0), mu_phi(4,0), met_phi(1,0);
131 
132  for (u_int i = 0; i<(*jets).size() && i<6; ++i) {
133  if ( parType_t( (*jets)[i].Et() ) <= p_minEt1 ) continue; //ET cut, leave NN inputs at default values (0)
134  jet_pt[i] = (*jets)[i].Et();
135  jet_eta[i] = (*jets)[i].eta();
136  jet_phi[i] = (*jets)[i].phi();
137  }
138  for (u_int i = 0; i < (*taus).size() && i<4; ++i) {
139  if ( parType_t( (*taus)[i].Et() ) <= p_minEt2 ) continue; //ET cut, leave NN inputs at default values (0)
140  tau_pt[i] = (*taus)[i].Et();
141  tau_eta[i] = (*taus)[i].eta();
142  tau_phi[i] = (*taus)[i].phi();
143  }
144  for (u_int i = 0; i < (*mus).size() && i<4; ++i) {
145  if ( parType_t( (*mus)[i].Et() ) <= p_minEt3 ) continue; //ET cut, leave NN inputs at default values (0)
146  mu_pt[i] = (*mus)[i].Et();
147  mu_eta[i] = (*mus)[i].eta();
148  mu_phi[i] = (*mus)[i].phi();
149  }
150  for (u_int i = 0; i < (*met).size() && i<1; ++i) {
151  if ( parType_t( (*met)[i].Et() ) <= p_minEt4 ) continue; //ET cut, leave NN inputs at default values (0)
152  met_pt[i] = (*met)[i].Et();
153  met_phi[i] = (*met)[i].phi();
154  }
155 
156 
157  TRG_MSG_DEBUG("Jet0: " << jet_pt[0] << ", " << jet_eta[0] << ", " << jet_phi[0] );
158  TRG_MSG_DEBUG("Jet1: " << jet_pt[1] << ", " << jet_eta[1] << ", " << jet_phi[1] );
159  TRG_MSG_DEBUG("Jet2: " << jet_pt[2] << ", " << jet_eta[2] << ", " << jet_phi[2] );
160  TRG_MSG_DEBUG("Jet3: " << jet_pt[3] << ", " << jet_eta[3] << ", " << jet_phi[3] );
161  TRG_MSG_DEBUG("Jet4: " << jet_pt[4] << ", " << jet_eta[4] << ", " << jet_phi[4] );
162  TRG_MSG_DEBUG("Jet5: " << jet_pt[5] << ", " << jet_eta[5] << ", " << jet_phi[5] );
163 
164  TRG_MSG_DEBUG("Tau0: " << tau_pt[0] << ", " << tau_eta[0] << ", " << tau_phi[0] );
165  TRG_MSG_DEBUG("Tau1: " << tau_pt[1] << ", " << tau_eta[1] << ", " << tau_phi[1] );
166  TRG_MSG_DEBUG("Tau2: " << tau_pt[2] << ", " << tau_eta[2] << ", " << tau_phi[2] );
167  TRG_MSG_DEBUG("Tau3: " << tau_pt[3] << ", " << tau_eta[3] << ", " << tau_phi[3] );
168 
169  TRG_MSG_DEBUG("Mu0: " << mu_pt[0] << ", " << mu_eta[0] << ", " << mu_phi[0] );
170  TRG_MSG_DEBUG("Mu1: " << mu_pt[1] << ", " << mu_eta[1] << ", " << mu_phi[1] );
171  TRG_MSG_DEBUG("Mu2: " << mu_pt[2] << ", " << mu_eta[2] << ", " << mu_phi[2] );
172  TRG_MSG_DEBUG("Mu3: " << mu_pt[3] << ", " << mu_eta[3] << ", " << mu_phi[3] );
173 
174  TRG_MSG_DEBUG("MET: " << met_pt[0] << ", " << met_phi[0] << std::endl);
175 
176  ADVAE2A::VAENetwork AD_Network( jet_pt[0], jet_eta[0], jet_phi[0],
177  jet_pt[1], jet_eta[1], jet_phi[1],
178  jet_pt[2], jet_eta[2], jet_phi[2],
179  jet_pt[3], jet_eta[3], jet_phi[3],
180  jet_pt[4], jet_eta[4], jet_phi[4],
181  jet_pt[5], jet_eta[5], jet_phi[5],
182  tau_pt[0], tau_eta[0], tau_phi[0],
183  tau_pt[1], tau_eta[1], tau_phi[1],
184  tau_pt[2], tau_eta[2], tau_phi[2],
185  tau_pt[3], tau_eta[3], tau_phi[3],
186  mu_pt [0], mu_eta [0], mu_phi [0],
187  mu_pt [1], mu_eta [1], mu_phi [1],
188  mu_pt [2], mu_eta [2], mu_phi [2],
189  mu_pt [3], mu_eta [3], mu_phi [3],
190  met_pt[0], met_phi[0] );
191  std::vector<int64_t> anomScoreInt64Vec = AD_Network.getAnomalyScoreInt64Vec();
192 
193  for(u_int i=0; i<numberOutputBits(); ++i) {
194  bool accept = false;
195  // Retrieve threshold
196  int32_t threshold = int32_t ( p_AnomalyScoreThresh[i] );
197  // Calculate event score
198  int64_t anomScoreInt64 = 0;
199  anomScoreInt64 = (p_ScaleSqr1[i] * anomScoreInt64Vec.at(0)*anomScoreInt64Vec.at(0) >> p_ScaleSqr_DropBits) +
200  (p_ScaleSqr2[i] * anomScoreInt64Vec.at(1)*anomScoreInt64Vec.at(1) >> p_ScaleSqr_DropBits) +
201  (p_ScaleSqr3[i] * anomScoreInt64Vec.at(2)*anomScoreInt64Vec.at(2) >> p_ScaleSqr_DropBits);
202  if ( anomScoreInt64 > threshold ) {
203  accept = true;
204  decision.setBit(i, true);
205  for ( u_int j = 0; j<6 && j<(*jets).size(); ++j ) output[i]->push_back((*jets)[j]);
206  for ( u_int j = 0; j<4 && j<(*taus).size(); ++j ) output[i]->push_back((*taus)[j]);
207  for ( u_int j = 0; j<4 && j<(*mus).size() ; ++j ) output[i]->push_back((*mus) [j]);
208  output[i]->push_back((*met)[0]);
209  }
210  output[i]->setAmbiguityFlag(hasAmbiguousInputs);
211 
212  if(fillHistos() and accept) {
213  fillHist1D(m_histAccept[i],anomScoreInt64);
214  } else if(fillHistos() && !accept) {
215  fillHist1D(m_histReject[i],anomScoreInt64);
216  }
217 
218  TRG_MSG_DEBUG("Decision for bit" << i << ": " << (accept?"pass":"fail") << " anomaly score = " << anomScoreInt64 << std::endl);
219  }
220  } else {
221  TCS_EXCEPTION("ADVAE_2A alg must have 4 inputs, but got " << input.size());
222  }
223 
224  return TCS::StatusCode::SUCCESS;
225 }
226 
227 TCS::StatusCode
228 TCS::ADVAE_2A::process( const std::vector<TCS::TOBArray const *> & input,
229  const std::vector<TCS::TOBArray *> & output,
230  Decision & decision )
231 {
232 
233 
234  if( input.size() == 4) {
235 
236  TCS::TOBArray const* jets = input[0];
237  TCS::TOBArray const* taus = input[1];
238  TCS::TOBArray const* mus = input[2];
239  TCS::TOBArray const* met = input[3];
240  TRG_MSG_DEBUG("Number of jets are " << (*jets).size());
241  TRG_MSG_DEBUG("Number of taus are " << (*taus).size());
242  TRG_MSG_DEBUG("Number of mus are " << (*mus).size());
243  TRG_MSG_DEBUG("Number of met are " << (*met).size());
244 
245  std::vector<u_int> jet_pt(6,0), tau_pt(4,0), mu_pt(4,0), met_pt(1,0);
246  std::vector<int> jet_eta(6,0), tau_eta(4,0), mu_eta(4,0); //no met_eta
247  std::vector<int> jet_phi(6,0), tau_phi(4,0), mu_phi(4,0), met_phi(1,0);
248 
249  for (u_int i = 0; i<(*jets).size() && i<6; ++i) {
250  if ( parType_t( (*jets)[i].Et() ) <= p_minEt1 ) continue; //ET cut, leave NN inputs at default values (0)
251  jet_pt[i] = (*jets)[i].Et();
252  jet_eta[i] = (*jets)[i].eta();
253  jet_phi[i] = (*jets)[i].phi();
254  }
255  for (u_int i = 0; i < (*taus).size() && i<4; ++i) {
256  if ( parType_t( (*taus)[i].Et() ) <= p_minEt2 ) continue; //ET cut, leave NN inputs at default values (0)
257  tau_pt[i] = (*taus)[i].Et();
258  tau_eta[i] = (*taus)[i].eta();
259  tau_phi[i] = (*taus)[i].phi();
260  }
261  for (u_int i = 0; i < (*mus).size() && i<4; ++i) {
262  if ( parType_t( (*mus)[i].Et() ) <= p_minEt3 ) continue; //ET cut, leave NN inputs at default values (0)
263  mu_pt[i] = (*mus)[i].Et();
264  mu_eta[i] = (*mus)[i].eta();
265  mu_phi[i] = (*mus)[i].phi();
266  }
267  for (u_int i = 0; i < (*met).size() && i<1; ++i) {
268  if ( parType_t( (*met)[i].Et() ) <= p_minEt4 ) continue; //ET cut, leave NN inputs at default values (0)
269  met_pt[i] = (*met)[i].Et();
270  met_phi[i] = (*met)[i].phi();
271  }
272 
273  ADVAE2A::VAENetwork AD_Network( jet_pt[0], jet_eta[0], jet_phi[0],
274  jet_pt[1], jet_eta[1], jet_phi[1],
275  jet_pt[2], jet_eta[2], jet_phi[2],
276  jet_pt[3], jet_eta[3], jet_phi[3],
277  jet_pt[4], jet_eta[4], jet_phi[4],
278  jet_pt[5], jet_eta[5], jet_phi[5],
279  tau_pt[0], tau_eta[0], tau_phi[0],
280  tau_pt[1], tau_eta[1], tau_phi[1],
281  tau_pt[2], tau_eta[2], tau_phi[2],
282  tau_pt[3], tau_eta[3], tau_phi[3],
283  mu_pt [0], mu_eta [0], mu_phi [0],
284  mu_pt [1], mu_eta [1], mu_phi [1],
285  mu_pt [2], mu_eta [2], mu_phi [2],
286  mu_pt [3], mu_eta [3], mu_phi [3],
287  met_pt[0], met_phi[0] );
288  std::vector<int64_t> anomScoreInt64Vec = AD_Network.getAnomalyScoreInt64Vec();
289 
290  for(u_int i=0; i<numberOutputBits(); ++i) {
291  bool accept = false;
292  // Retrieve threshold
293  int32_t threshold = int32_t ( p_AnomalyScoreThresh[i] );
294  // Calculate event score
295  int64_t anomScoreInt64 = 0;
296  anomScoreInt64 = p_ScaleSqr1[i]/std::pow(2,p_ScaleSqr_DropBits) * anomScoreInt64Vec.at(0)*anomScoreInt64Vec.at(0) +
297  p_ScaleSqr2[i]/std::pow(2,p_ScaleSqr_DropBits) * anomScoreInt64Vec.at(1)*anomScoreInt64Vec.at(1) +
298  p_ScaleSqr3[i]/std::pow(2,p_ScaleSqr_DropBits) * anomScoreInt64Vec.at(2)*anomScoreInt64Vec.at(2);
299  if ( anomScoreInt64 > threshold ) {
300  accept = true;
301  decision.setBit(i, true);
302  for ( u_int j = 0; j<6 && j<(*jets).size(); ++j ) output[i]->push_back((*jets)[j]);
303  for ( u_int j = 0; j<4 && j<(*taus).size(); ++j ) output[i]->push_back((*taus)[j]);
304  for ( u_int j = 0; j<4 && j<(*mus).size() ; ++j ) output[i]->push_back((*mus) [j]);
305  output[i]->push_back((*met)[0]);
306  }
307 
308  if(fillHistos() and accept) {
309  fillHist1D(m_histAccept[i],anomScoreInt64);
310  } else if(fillHistos() && !accept) {
311  fillHist1D(m_histReject[i],anomScoreInt64);
312  }
313 
314  TRG_MSG_DEBUG("Decision for bit" << i << ": " << (accept?"pass":"fail") << " anomaly score = " << anomScoreInt64 << std::endl);
315  }
316  } else {
317  TCS_EXCEPTION("ADVAE_2A alg must have 4 inputs, but got " << input.size());
318  }
319 
320  return TCS::StatusCode::SUCCESS;
321 }
TCS::ADVAE_2A::initialize
virtual StatusCode initialize()
Definition: AnomDetVAE.cxx:66
TCS::StatusCode::SUCCESS
@ SUCCESS
Definition: Trigger/TrigT1/L1Topo/L1TopoCommon/L1TopoCommon/StatusCode.h:17
TCS::parType_t
uint32_t parType_t
Definition: Parameter.h:22
CutsMETMaker::accept
StatusCode accept(const xAOD::Muon *mu)
Definition: CutsMETMaker.cxx:18
defineDB.jets
jets
Definition: JetTagCalibration/share/defineDB.py:24
athena.value
value
Definition: athena.py:124
TCS::ADVAE_2A
Definition: AnomDetVAE.h:16
Decision.h
const
bool const RAWDATA *ch2 const
Definition: LArRodBlockPhysicsV0.cxx:560
TCS::DecisionAlg
Definition: Trigger/TrigT1/L1Topo/L1TopoInterfaces/L1TopoInterfaces/DecisionAlg.h:22
TCS::Decision::setBit
void setBit(unsigned int index, bool value)
Definition: L1Topo/L1TopoInterfaces/Root/Decision.cxx:12
met
Definition: IMETSignificance.h:24
lumiFormat.i
int i
Definition: lumiFormat.py:85
TCS_EXCEPTION
#define TCS_EXCEPTION(MSG)
Definition: Trigger/TrigT1/L1Topo/L1TopoCommon/L1TopoCommon/Exception.h:14
TCS::TOBArray
Definition: TOBArray.h:24
plotIsoValidation.mu_phi
mu_phi
Definition: plotIsoValidation.py:152
TCS::Decision
Definition: L1Topo/L1TopoInterfaces/L1TopoInterfaces/Decision.h:19
TRG_MSG_INFO
#define TRG_MSG_INFO(x)
Definition: Trigger/TrigConfiguration/TrigConfBase/TrigConfBase/MsgStreamMacros.h:27
TCS::ADVAE_2A::processBitCorrect
virtual StatusCode processBitCorrect(const std::vector< TCS::TOBArray const * > &input, const std::vector< TCS::TOBArray * > &output, Decision &decison)
Definition: AnomDetVAE.cxx:105
TCS::ADVAE_2A::process
virtual StatusCode process(const std::vector< TCS::TOBArray const * > &input, const std::vector< TCS::TOBArray * > &output, Decision &decison)
Definition: AnomDetVAE.cxx:228
REGISTER_ALG_TCS
#define REGISTER_ALG_TCS(CLASS)
Definition: AlgFactory.h:62
AnomDetVAE.h
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
ActsTrk::to_string
std::string to_string(const DetectorType &type)
Definition: GeometryDefs.h:34
threshold
Definition: chainparser.cxx:74
TCS
Definition: Global/GlobalSimulation/src/IO/Decision.h:18
TSU::isAmbiguousAnywhere
bool isAmbiguousAnywhere(TCS::TOBArray const *tobs, size_t pos, unsigned minEt=0)
Definition: Trigger/TrigT1/L1Topo/L1TopoSimulationUtils/Root/Helpers.cxx:37
plotIsoValidation.mu_eta
mu_eta
Definition: plotIsoValidation.py:151
Trk::jet_phi
@ jet_phi
Definition: JetVtxParamDefs.h:28
Helpers.h
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
Exception.h
TCS::ADVAE_2A::~ADVAE_2A
virtual ~ADVAE_2A()
Definition: AnomDetVAE.cxx:62
TRG_MSG_DEBUG
#define TRG_MSG_DEBUG(x)
Definition: Trigger/TrigConfiguration/TrigConfBase/TrigConfBase/MsgStreamMacros.h:25
TCS::StatusCode
Definition: Trigger/TrigT1/L1Topo/L1TopoCommon/L1TopoCommon/StatusCode.h:15
plotIsoValidation.mu_pt
mu_pt
Definition: plotIsoValidation.py:150
Generated on Mon Sep 29 2025 21:06:30 for ATLAS Offline Software by   doxygen 1.8.18