ATLAS Offline Software
TFCSGANEtaSlice.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // TFCSGANEtaSlice.cxx, (c) ATLAS Detector software //
8 
9 // class header include
11 
13 
14 #include "CLHEP/Random/RandGauss.h"
15 
16 #include "TFitResult.h"
17 #include "TFile.h"
18 #include "TTree.h"
19 #include "TH1D.h"
20 #include "TMath.h"
21 
22 #include <iostream>
23 #include <fstream>
24 #include <string>
25 #include <sstream>
26 #include <cmath>
27 
28 
30 
32  const TFCSGANXMLParameters &param)
33  : m_pid(pid), m_etaMin(etaMin), m_etaMax(etaMax), m_param(param) {}
34 
36  // Deleting a nullptr is a noop
37  delete m_gan_all;
38  delete m_gan_low;
39  delete m_gan_high;
40 }
41 
43  if (m_net_all != nullptr)
44  return m_net_all.get();
45  return m_gan_all;
46 }
48  if (m_net_low != nullptr)
49  return m_net_low.get();
50  return m_gan_low;
51 }
53  if (m_net_high != nullptr)
54  return m_net_high.get();
55  return m_gan_high;
56 }
57 
59  if (m_pid == 211 || m_pid == 2212) {
60  if (GetNetAll() == nullptr) {
61  return false;
62  }
63  } else {
64  if (GetNetHigh() == nullptr || GetNetLow() == nullptr) {
65  return false;
66  }
67  }
68  return true;
69 }
70 
72  // Now load new data
73  std::string inputFileName;
74 
77 
78  bool success = true;
79 
80  if (m_pid == 211) {
81  inputFileName = m_param.GetInputFolder() + "/neural_net_" +
83  "_" + std::to_string(m_etaMax) + "_All.*";
84  ATH_MSG_DEBUG("Gan input file name " << inputFileName);
85  m_net_all = TFCSNetworkFactory::create(inputFileName);
86  if (m_net_all == nullptr)
87  success = false;
88  } else if (m_pid == 2212) {
89  inputFileName = m_param.GetInputFolder() + "/neural_net_" +
91  "_" + std::to_string(m_etaMax) + "_High10.*";
92  ATH_MSG_DEBUG("Gan input file name " << inputFileName);
93  m_net_all = TFCSNetworkFactory::create(inputFileName);
94  if (m_net_all == nullptr)
95  success = false;
96  } else {
97  inputFileName = m_param.GetInputFolder() + "/neural_net_" +
99  "_" + std::to_string(m_etaMax) + "_High12.*";
100  ATH_MSG_DEBUG("Gan input file name " << inputFileName);
101  m_net_high = TFCSNetworkFactory::create(inputFileName);
102  if (m_net_high == nullptr)
103  success = false;
104 
105  inputFileName = m_param.GetInputFolder() + "/neural_net_" +
107  "_" + std::to_string(m_etaMax) + "_UltraLow12.*";
108  m_net_low = TFCSNetworkFactory::create(inputFileName);
109  if (m_net_low == nullptr)
110  success = false;
111  }
112  return success;
113 }
114 
116  std::string rootFileName = m_param.GetInputFolder() + "/rootFiles/pid" +
117  std::to_string(m_pid) + "_E1048576_eta_" +
118  std::to_string(m_etaMin) + "_" +
119  std::to_string(m_etaMin + 5) + ".root";
120  ATH_MSG_DEBUG("Opening file " << rootFileName);
121  TFile *file = TFile::Open(rootFileName.c_str(), "read");
122  for (int layer : m_param.GetRelevantLayers()) {
123  ATH_MSG_DEBUG("Layer " << layer);
125  TH2D *h2 = &binsInLayers[layer];
126 
127  std::string histoName = "r" + std::to_string(layer) + "w";
128  TH1D *h1 = (TH1D *)file->Get(histoName.c_str());
129  if (TMath::IsNaN(h1->Integral())) {
130  histoName = "r" + std::to_string(layer);
131  h1 = (TH1D *)file->Get(histoName.c_str());
132  }
133 
134  TAxis *x = (TAxis *)h2->GetXaxis();
135  for (int ix = 1; ix <= h2->GetNbinsX(); ++ix) {
136  ATH_MSG_DEBUG(ix);
137  h1->GetXaxis()->SetRangeUser(x->GetBinLowEdge(ix), x->GetBinUpEdge(ix));
138 
139  double result = 0;
140  if (h1->Integral() > 0 && h1->GetNbinsX() > 2) {
141  TFitResultPtr res(0);
142 
143  res = h1->Fit("expo", "SQ");
144  if (res >= 0 && !std::isnan(res->Parameter(0))) {
145  result = res->Parameter(1);
146  }
147  }
148  m_allFitResults[layer].push_back(result);
149  }
150  }
151  ATH_MSG_DEBUG("Done initialisaing fits");
152 }
153 
155  std::string rootFileName = m_param.GetInputFolder() + "/rootFiles/pid" +
156  std::to_string(m_pid) + "_E65536_eta_" +
157  std::to_string(m_etaMin) + "_" +
158  std::to_string(m_etaMin + 5) + "_validation.root";
159  ATH_MSG_DEBUG("Opening file " << rootFileName);
160  TFile *file = TFile::Open(rootFileName.c_str(), "read");
161  for (int layer : m_param.GetRelevantLayers()) {
162  std::string branchName = "extrapWeight_" + std::to_string(layer);
163  TH1D *h = new TH1D("h", "h", 100, 0.01, 1);
164  TTree *tree = (TTree *)file->Get("rootTree");
165  std::string command = branchName + ">>h";
166  tree->Draw(command.c_str());
167  m_extrapolatorWeights[layer] = h->GetMean();
168  ATH_MSG_DEBUG("Extrapolation: layer " << layer << " mean "
170  }
171 }
172 
176  TFCSSimulationState simulstate) const {
177  double randUniformZ = 0.;
179 
180  int maxExp = 0, minExp = 0;
181  if (m_pid == 22 || std::abs(m_pid) == 11) {
182  if (truth->P() >
183  4096) { // This is the momentum, not the energy, because the split is
184  // based on the samples which are produced with the momentum
185  maxExp = 22;
186  minExp = 12;
187  } else {
188  maxExp = 12;
189  minExp = 6;
190  }
191  } else if (std::abs(m_pid) == 211) {
192  maxExp = 22;
193  minExp = 8;
194  } else if (std::abs(m_pid) == 2212) {
195  maxExp = 22;
196  minExp = 10;
197  }
198 
199  int p_min = std::pow(2, minExp);
200  int p_max = std::pow(2, maxExp);
201  // Keep min and max without mass offset as we do not train on antiparticles
202  double Ekin_min =
203  std::sqrt(std::pow(p_min, 2) + std::pow(truth->M(), 2)) - truth->M();
204  double Ekin_max =
205  std::sqrt(std::pow(p_max, 2) + std::pow(truth->M(), 2)) - truth->M();
206 
207  for (int i = 0; i < m_param.GetLatentSpaceSize(); i++) {
208  randUniformZ = CLHEP::RandGauss::shoot(simulstate.randomEngine(), 0.5, 0.5);
209  inputs["Noise"].insert(std::pair<std::string, double>(
210  "variable_" + std::to_string(i), randUniformZ));
211  }
212 
213  // double e = log(truth->Ekin()/Ekin_min)/log(Ekin_max/Ekin_min) ;
214  // Could be uncommented , but would need the line above too
215  // ATH_MSG_DEBUG( "Check label: " << e <<" Ekin:" << truth->Ekin() <<" p:" <<
216  // truth->P() <<" mass:" << truth->M() <<" Ekin_off:" <<
217  // truth->Ekin_off() << " Ekin_min:"<<Ekin_min<<"
218  // Ekin_max:"<<Ekin_max);
219  // inputs["mycond"].insert ( std::pair<std::string,double>("variable_0",
220  // truth->Ekin()/(std::pow(2,maxExp))) ); //Old conditioning using linear
221  // interpolation, now use logaritminc interpolation
222  inputs["mycond"].insert(std::pair<std::string, double>(
223  "variable_0", log(truth->Ekin() / Ekin_min) / log(Ekin_max / Ekin_min)));
224 
225  if (m_param.GetGANVersion() >= 2) {
226  if (false) { // conditioning on eta, should only be needed in transition
227  // regions and added only to the GANs that use it, for now all
228  // GANs have 3 conditioning inputs so filling zeros
229  inputs["mycond"].insert(std::pair<std::string, double>(
230  "variable_1", std::abs(extrapol->IDCaloBoundary_eta())));
231  } else {
232  inputs["mycond"].insert(std::pair<std::string, double>("variable_1", 0));
233  }
234  }
235 
237  if (m_param.GetGANVersion() == 1 || m_pid == 211 || m_pid == 2212) {
239  } else {
240  if (truth->P() >
241  4096) { // This is the momentum, not the energy, because the split is
242  // based on the samples which are produced with the momentum
243  ATH_MSG_DEBUG("Computing outputs given inputs for high");
245  } else {
247  }
248  }
249  ATH_MSG_DEBUG("Start Network inputs ~~~~~~~~");
251  ATH_MSG_DEBUG("End Network inputs ~~~~~~~~");
252  ATH_MSG_DEBUG("Start Network outputs ~~~~~~~~");
254  ATH_MSG_DEBUG("End Network outputs ~~~~~~~~");
255  return outputs;
256 }
257 
259  ATH_MSG_INFO("LWTNN Handler parameters");
260  ATH_MSG_INFO(" pid: " << m_pid);
261  ATH_MSG_INFO(" etaMin:" << m_etaMin);
262  ATH_MSG_INFO(" etaMax: " << m_etaMax);
263  m_param.Print();
264 }
TFCSGANXMLParameters::Print
void Print() const
Definition: TFCSGANXMLParameters.cxx:124
VNetworkBase::NetworkOutputs
std::map< std::string, double > NetworkOutputs
Format for network outputs.
Definition: VNetworkBase.h:100
TFCSGANEtaSlice::LoadGAN
bool LoadGAN()
Definition: TFCSGANEtaSlice.cxx:71
get_generator_info.result
result
Definition: get_generator_info.py:21
TFCSGANEtaSlice::m_allFitResults
FitResultsPerLayer m_allFitResults
Definition: TFCSGANEtaSlice.h:63
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
TFCSGANEtaSlice::m_net_low
std::unique_ptr< VNetworkBase > m_net_low
Definition: TFCSGANEtaSlice.h:73
VNetworkBase::representNetworkOutputs
static std::string representNetworkOutputs(NetworkOutputs const &outputs, int maxValues=3)
String representation of network outputs.
Definition: VNetworkBase.cxx:57
tree
TChain * tree
Definition: tile_monitor.h:30
TFCSGANEtaSlice::m_net_all
std::unique_ptr< VNetworkBase > m_net_all
Definition: TFCSGANEtaSlice.h:72
xAOD::etaMax
etaMax
Definition: HIEventShape_v2.cxx:46
TFCSGANXMLParameters::Binning
std::map< int, TH2D > Binning
Definition: TFCSGANXMLParameters.h:23
TFCSGANEtaSlice::m_pid
int m_pid
Definition: TFCSGANEtaSlice.h:57
TFCSGANEtaSlice::m_etaMin
int m_etaMin
Definition: TFCSGANEtaSlice.h:58
TFCSGANXMLParameters::GetRelevantLayers
const std::vector< int > & GetRelevantLayers() const
Definition: TFCSGANXMLParameters.h:32
TFCSExtrapolationState
Definition: TFCSExtrapolationState.h:13
TFCSSimulationState::randomEngine
CLHEP::HepRandomEngine * randomEngine()
Definition: TFCSSimulationState.h:36
RunActsMaterialValidation.extrapol
extrapol
Definition: RunActsMaterialValidation.py:91
read_hist_ntuple.h1
h1
Definition: read_hist_ntuple.py:21
TFCSGANXMLParameters
Definition: TFCSGANXMLParameters.h:21
VNetworkBase::compute
virtual NetworkOutputs compute(NetworkInputs const &inputs) const =0
Function to pass values to the network.
x
#define x
TFCSNetworkFactory::create
static std::unique_ptr< VNetworkBase > create(std::string input)
Given a string, make a network.
Definition: TFCSNetworkFactory.cxx:66
postInclude.inputs
inputs
Definition: postInclude.SortInput.py:15
TFCSGANEtaSlice::GetNetAll
VNetworkBase * GetNetAll() const
Definition: TFCSGANEtaSlice.cxx:42
lumiFormat.i
int i
Definition: lumiFormat.py:85
TFCSGANEtaSlice::GetNetHigh
VNetworkBase * GetNetHigh() const
Definition: TFCSGANEtaSlice.cxx:52
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
TFCSGANEtaSlice::m_param
TFCSGANXMLParameters m_param
Definition: TFCSGANEtaSlice.h:83
TFCSTruthState::Ekin
double Ekin() const
Definition: TFCSTruthState.h:26
TRT::Hit::layer
@ layer
Definition: HitInfo.h:79
res
std::pair< std::vector< unsigned int >, bool > res
Definition: JetGroupProductTest.cxx:14
ParticleGun_EoverP_Config.pid
pid
Definition: ParticleGun_EoverP_Config.py:62
file
TFile * file
Definition: tile_monitor.h:29
TFCSGANEtaSlice::TFCSGANEtaSlice
TFCSGANEtaSlice()
Definition: TFCSGANEtaSlice.cxx:29
TFCSGANEtaSlice::Print
void Print() const
Definition: TFCSGANEtaSlice.cxx:258
TFCSGANEtaSlice::CalculateMeanPointFromDistributionOfR
void CalculateMeanPointFromDistributionOfR()
Definition: TFCSGANEtaSlice.cxx:115
TFCSGANEtaSlice.h
TFCSGANEtaSlice::m_gan_low
TFCSGANLWTNNHandler * m_gan_low
Definition: TFCSGANEtaSlice.h:69
TFCSGANEtaSlice::IsGanCorrectlyLoaded
bool IsGanCorrectlyLoaded() const
Definition: TFCSGANEtaSlice.cxx:58
python.CreateTierZeroArgdict.outputs
outputs
Definition: CreateTierZeroArgdict.py:189
TFCSGANEtaSlice::GetNetLow
VNetworkBase * GetNetLow() const
Definition: TFCSGANEtaSlice.cxx:47
TFCSGANEtaSlice::m_extrapolatorWeights
ExtrapolatorWeights m_extrapolatorWeights
Definition: TFCSGANEtaSlice.h:64
TFCSGANEtaSlice::NetworkInputs
std::map< std::string, std::map< std::string, double > > NetworkInputs
Definition: TFCSGANEtaSlice.h:39
VNetworkBase
Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration.
Definition: VNetworkBase.h:38
TFCSGANEtaSlice::m_net_high
std::unique_ptr< VNetworkBase > m_net_high
Definition: TFCSGANEtaSlice.h:74
ActsTrk::to_string
std::string to_string(const DetectorType &type)
Definition: GeometryDefs.h:34
TFCSGANXMLParameters::GetLatentSpaceSize
int GetLatentSpaceSize() const
Definition: TFCSGANXMLParameters.h:34
VNetworkBase::representNetworkInputs
static std::string representNetworkInputs(NetworkInputs const &inputs, int maxValues=3)
String representation of network inputs.
Definition: VNetworkBase.cxx:37
TFCSGANXMLParameters::GetGANVersion
int GetGANVersion() const
Definition: TFCSGANXMLParameters.h:35
TFCSGANEtaSlice::m_etaMax
int m_etaMax
Definition: TFCSGANEtaSlice.h:59
TFCSGANEtaSlice::ExtractExtrapolatorMeansFromInputs
void ExtractExtrapolatorMeansFromInputs()
Definition: TFCSGANEtaSlice.cxx:154
TFCSGANEtaSlice::GetNetworkOutputs
NetworkOutputs GetNetworkOutputs(const TFCSTruthState *truth, const TFCSExtrapolationState *extrapol, TFCSSimulationState simulstate) const
Definition: TFCSGANEtaSlice.cxx:174
TFCSGANXMLParameters::GetInputFolder
const std::string & GetInputFolder() const
Definition: TFCSGANXMLParameters.h:37
LArCellBinning.etaMin
etaMin
Definition: LArCellBinning.py:84
h
TFCSGANEtaSlice::~TFCSGANEtaSlice
virtual ~TFCSGANEtaSlice()
Definition: TFCSGANEtaSlice.cxx:35
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
TFCSGANEtaSlice::m_gan_high
TFCSGANLWTNNHandler * m_gan_high
Definition: TFCSGANEtaSlice.h:70
TFCSNetworkFactory.h
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
TFCSTruthState
Definition: TFCSTruthState.h:13
TFCSSimulationState
Definition: TFCSSimulationState.h:32
get_generator_info.command
string command
Definition: get_generator_info.py:38
TFCSGANEtaSlice::m_gan_all
TFCSGANLWTNNHandler * m_gan_all
Definition: TFCSGANEtaSlice.h:68
TFCSGANXMLParameters::GetBinning
const Binning & GetBinning() const
Definition: TFCSGANXMLParameters.h:33