ATLAS Offline Software
AnomDetVAE.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 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 AnomalyScoreThresh
13 **********************************/
14 
15 #include <cmath>
16 
17 #include "L1TopoAlgorithms/AnomDetVAE.h"
18 #include "L1TopoCommon/Exception.h"
19 #include "L1TopoInterfaces/Decision.h"
20 
21 REGISTER_ALG_TCS(ADVAE_2A)
22 
23 
24 TCS::ADVAE_2A::ADVAE_2A(const std::string & name) : DecisionAlg(name)
25 {
26  defineParameter("InputWidth1", 6);
27  defineParameter("InputWidth2", 6);
28  defineParameter("InputWidth3", 6);
29  defineParameter("InputWidth4", 1);
30  defineParameter("MaxTob1", 6);
31  defineParameter("MaxTob2", 4);
32  defineParameter("MaxTob3", 4);
33  defineParameter("MaxTob4", 1);
34  defineParameter("NumResultBits", 2);
35 
36  // Version parameter, used for L1TopoFW book-keeping, no practical application
37  defineParameter("ADVAEVersion", 1);
38  //The value used is AD threshold * 3 (since there are 3 gaussians) * 1024 (to make the decimal points available)
39  defineParameter("AnomalyScoreThresh", 3875, 0);
40  defineParameter("AnomalyScoreThresh", 3875, 1);
41 
42  setNumberOutputBits(2);
43 }
44 
45 TCS::ADVAE_2A::~ADVAE_2A(){}
46 
47 
48 TCS::StatusCode
49 TCS::ADVAE_2A::initialize() {
50  p_NumberLeading1 = parameter("InputWidth1").value();
51  p_NumberLeading2 = parameter("InputWidth2").value();
52  p_NumberLeading3 = parameter("InputWidth3").value();
53  p_NumberLeading4 = parameter("InputWidth4").value();
54 
55  if(parameter("MaxTob1").value() > 0) p_NumberLeading1 = parameter("MaxTob1").value();
56  if(parameter("MaxTob2").value() > 0) p_NumberLeading2 = parameter("MaxTob2").value();
57  if(parameter("MaxTob3").value() > 0) p_NumberLeading3 = parameter("MaxTob3").value();
58  if(parameter("MaxTob4").value() > 0) p_NumberLeading4 = parameter("MaxTob4").value();
59 
60  for(unsigned int i=0; i<numberOutputBits(); ++i) {
61  p_AnomalyScoreThresh[i] = parameter("AnomalyScoreThresh", i).value();
62  }
63 
64  TRG_MSG_INFO("number output : " << numberOutputBits());
65 
66  // book histograms
67  for(unsigned int i=0; i<numberOutputBits(); ++i) {
68  std::string hname_accept = "hAnomalyScore_accept_bit"+std::to_string((int)i);
69  std::string hname_reject = "hAnomalyScore_reject_bit"+std::to_string((int)i);
70  // score
71  bookHist(m_histAccept, hname_accept, "ADScore", 500, 0, 5000);
72  bookHist(m_histReject, hname_reject, "ADScore", 500, 0, 5000);
73  }
74 
75  return StatusCode::SUCCESS;
76 }
77 
78 
79 
80 TCS::StatusCode
81 TCS::ADVAE_2A::processBitCorrect( const std::vector<TCS::TOBArray const *> & input,
82  const std::vector<TCS::TOBArray *> & output,
83  Decision & decision )
84 {
85 
86 
87  if( input.size() == 4) {
88 
89  TCS::TOBArray const* jets = input[0];
90  TCS::TOBArray const* taus = input[1];
91  TCS::TOBArray const* mus = input[2];
92  TCS::TOBArray const* met = input[3];
93  TRG_MSG_DEBUG("Number of jets are " << (*jets).size());
94  TRG_MSG_DEBUG("Number of taus are " << (*taus).size());
95  TRG_MSG_DEBUG("Number of mus are " << (*mus).size());
96  TRG_MSG_DEBUG("Number of met are " << (*met).size());
97 
98  std::vector<u_int> jet_pt(6,0), tau_pt(4,0), mu_pt(4,0);
99  std::vector<int> jet_eta(6,0), tau_eta(4,0), mu_eta(4,0);
100  std::vector<int> jet_phi(6,0), tau_phi(4,0), mu_phi(4,0);
101 
102  for (u_int i = 0; i<(*jets).size() && i<6; ++i) {
103  jet_pt[i] = (*jets)[i].Et();
104  jet_eta[i] = (*jets)[i].eta();
105  jet_phi[i] = (*jets)[i].phi();
106  }
107  for (u_int i = 0; i < (*taus).size() && i<4; ++i) {
108  tau_pt[i] = (*taus)[i].Et();
109  tau_eta[i] = (*taus)[i].eta();
110  tau_phi[i] = (*taus)[i].phi();
111  }
112  for (u_int i = 0; i < (*mus).size() && i<4; ++i) {
113  mu_pt[i] = (*mus)[i].Et();
114  mu_eta[i] = (*mus)[i].eta();
115  mu_phi[i] = (*mus)[i].phi();
116  }
117 
119  /*
120  Trig::ADScore AD_Score( jet_pt[0], jet_eta[0], jet_phi[0],
121  jet_pt[1], jet_eta[1], jet_phi[1],
122  jet_pt[2], jet_eta[2], jet_phi[2],
123  jet_pt[3], jet_eta[3], jet_phi[3],
124  jet_pt[4], jet_eta[4], jet_phi[4],
125  jet_pt[5], jet_eta[5], jet_phi[5],
126  tau_pt[0], tau_eta[0], tau_phi[0],
127  tau_pt[1], tau_eta[1], tau_phi[1],
128  tau_pt[2], tau_eta[2], tau_phi[2],
129  tau_pt[3], tau_eta[3], tau_phi[3],
130  mu_pt [0], mu_eta [0], mu_phi [0],
131  mu_pt [1], mu_eta [1], mu_phi [1],
132  mu_pt [2], mu_eta [2], mu_phi [2],
133  mu_pt [3], mu_eta [3], mu_phi [3],
134  (*met)[0].Et(), (*met)[0].phi() );
135  auto scoreVec = AD_Score.myScore();
136  */
137  std::vector<double> scoreVec(3,0);
138 
139  TRG_MSG_DEBUG("The anomaly score is " << scoreVec[0] << ", " << scoreVec[1] << ", " << scoreVec[2] << std::endl);
140  parType_t score = parType_t ( (scoreVec[0]*scoreVec[0] + scoreVec[1]*scoreVec[1] + scoreVec[2]*scoreVec[2])*1024 );
141 
142  for(u_int i=0; i<numberOutputBits(); ++i) {
143  bool accept = false;
144  if ( score > p_AnomalyScoreThresh[i] ) {
145  accept = true;
146  decision.setBit(i, true);
147  for ( u_int j = 0; j<6 && j<(*jets).size(); ++j ) output[i]->push_back((*jets)[j]);
148  for ( u_int j = 0; j<4 && j<(*taus).size(); ++j ) output[i]->push_back((*taus)[j]);
149  for ( u_int j = 0; j<4 && j<(*mus).size() ; ++j ) output[i]->push_back((*mus) [j]);
150  output[i]->push_back((*met)[0]);
151  }
152 
153  if(fillHistos() and accept) {
154  fillHist1D(m_histAccept[i],score);
155  } else if(fillHistos() && !accept) {
156  fillHist1D(m_histReject[i],score);
157  }
158 
159  TRG_MSG_DEBUG("Decision for bit" << i << ": " << (accept?"pass":"fail") << " anomaly score = " << score << std::endl);
160  }
161  } else {
162  TCS_EXCEPTION("ADVAE_2A alg must have 4 inputs, but got " << input.size());
163  }
164 
165  return TCS::StatusCode::SUCCESS;
166 }
167 
168 TCS::StatusCode
169 TCS::ADVAE_2A::process( const std::vector<TCS::TOBArray const *> & input,
170  const std::vector<TCS::TOBArray *> & output,
171  Decision & decision )
172 {
173 
174 
175  if( input.size() == 4) {
176 
177  TCS::TOBArray const* jets = input[0];
178  TCS::TOBArray const* taus = input[1];
179  TCS::TOBArray const* mus = input[2];
180  TCS::TOBArray const* met = input[3];
181  TRG_MSG_DEBUG("Number of jets are " << (*jets).size());
182  TRG_MSG_DEBUG("Number of taus are " << (*taus).size());
183  TRG_MSG_DEBUG("Number of mus are " << (*mus).size());
184  TRG_MSG_DEBUG("Number of met are " << (*met).size());
185 
186  std::vector<u_int> jet_pt(6,0), tau_pt(4,0), mu_pt(4,0);
187  std::vector<int> jet_eta(6,0), tau_eta(4,0), mu_eta(4,0);
188  std::vector<int> jet_phi(6,0), tau_phi(4,0), mu_phi(4,0);
189 
190  for (u_int i = 0; i<(*jets).size() && i<6; ++i) {
191  jet_pt[i] = (*jets)[i].Et();
192  jet_eta[i] = (*jets)[i].eta();
193  jet_phi[i] = (*jets)[i].phi();
194  }
195  for (u_int i = 0; i < (*taus).size() && i<4; ++i) {
196  tau_pt[i] = (*taus)[i].Et();
197  tau_eta[i] = (*taus)[i].eta();
198  tau_phi[i] = (*taus)[i].phi();
199  }
200  for (u_int i = 0; i < (*mus).size() && i<4; ++i) {
201  mu_pt[i] = (*mus)[i].Et();
202  mu_eta[i] = (*mus)[i].eta();
203  mu_phi[i] = (*mus)[i].phi();
204  }
206  /*
207  Trig::ADScore AD_Score( jet_pt[0], jet_eta[0], jet_phi[0],
208  jet_pt[1], jet_eta[1], jet_phi[1],
209  jet_pt[2], jet_eta[2], jet_phi[2],
210  jet_pt[3], jet_eta[3], jet_phi[3],
211  jet_pt[4], jet_eta[4], jet_phi[4],
212  jet_pt[5], jet_eta[5], jet_phi[5],
213  tau_pt[0], tau_eta[0], tau_phi[0],
214  tau_pt[1], tau_eta[1], tau_phi[1],
215  tau_pt[2], tau_eta[2], tau_phi[2],
216  tau_pt[3], tau_eta[3], tau_phi[3],
217  mu_pt [0], mu_eta [0], mu_phi [0],
218  mu_pt [1], mu_eta [1], mu_phi [1],
219  mu_pt [2], mu_eta [2], mu_phi [2],
220  mu_pt [3], mu_eta [3], mu_phi [3],
221  (*met)[0].Et(), (*met)[0].phi() );
222  auto scoreVec = AD_Score.myScore();
223  */
224  std::vector<double> scoreVec(3,0);
225 
226  TRG_MSG_DEBUG("The anomaly score is " << scoreVec[0] << ", " << scoreVec[1] << ", " << scoreVec[2] << std::endl);
227  parType_t score = parType_t ( (scoreVec[0]*scoreVec[0] + scoreVec[1]*scoreVec[1] + scoreVec[2]*scoreVec[2])*1024 );
228 
229  for(u_int i=0; i<numberOutputBits(); ++i) {
230  bool accept = false;
231  if ( score > p_AnomalyScoreThresh[i] ) {
232  accept = true;
233  decision.setBit(i, true);
234  for ( u_int j = 0; j<6 && j<(*jets).size(); ++j ) output[i]->push_back((*jets)[j]);
235  for ( u_int j = 0; j<4 && j<(*taus).size(); ++j ) output[i]->push_back((*taus)[j]);
236  for ( u_int j = 0; j<4 && j<(*mus).size() ; ++j ) output[i]->push_back((*mus) [j]);
237  output[i]->push_back((*met)[0]);
238  }
239 
240  if(fillHistos() and accept) {
241  fillHist1D(m_histAccept[i],score);
242  } else if(fillHistos() && !accept) {
243  fillHist1D(m_histReject[i],score);
244  }
245 
246  TRG_MSG_DEBUG("Decision for bit" << i << ": " << (accept?"pass":"fail") << " anomaly score = " << score << std::endl);
247  }
248  } else {
249  TCS_EXCEPTION("ADVAE_2A alg must have 4 inputs, but got " << input.size());
250  }
251 
252  return TCS::StatusCode::SUCCESS;
253 }
TCS::ADVAE_2A::initialize
virtual StatusCode initialize()
Definition: AnomDetVAE.cxx:49
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
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
PlotPulseshapeFromCool.input
input
Definition: PlotPulseshapeFromCool.py:106
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
merge.output
output
Definition: merge.py:17
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:81
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:169
REGISTER_ALG_TCS
#define REGISTER_ALG_TCS(CLASS)
Definition: AlgFactory.h:62
AnomDetVAE.h
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
ActsTrk::to_string
std::string to_string(const DetectorType &type)
Definition: GeometryDefs.h:34
xAOD::score
@ score
Definition: TrackingPrimitives.h:513
TCS
Definition: Global/GlobalSimulation/src/IO/Decision.h:18
plotIsoValidation.mu_eta
mu_eta
Definition: plotIsoValidation.py:151
defineDB.jets
list jets
Definition: JetTagCalibration/share/defineDB.py:24
Trk::jet_phi
@ jet_phi
Definition: JetVtxParamDefs.h:28
Exception.h
TCS::ADVAE_2A::~ADVAE_2A
virtual ~ADVAE_2A()
Definition: AnomDetVAE.cxx:45
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 Sun Dec 22 2024 21:06:56 for ATLAS Offline Software by   doxygen 1.8.18