ATLAS Offline Software
TFCSEnergyAndHitGANV2.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 
7 #include "CLHEP/Random/RandFlat.h"
8 #include "CLHEP/Random/RandGauss.h"
9 #include "HepPDT/ParticleData.hh"
10 #include "HepPDT/ParticleDataTable.hh"
17 #include "TF1.h"
18 #include "TFile.h"
19 #include "TH2D.h"
20 
21 #if defined(__FastCaloSimStandAlone__)
22 #include "CLHEP/Random/TRandomEngine.h"
23 #else
24 #include <CLHEP/Random/RanluxEngine.h>
25 #endif
26 
27 #include <fstream>
28 #include <iostream>
29 #include <limits>
30 
31 //=============================================
32 //======= TFCSEnergyAndHitGANV2 =========
33 //=============================================
34 
36  const char *title)
39 }
40 
42  if (m_slice != nullptr) {
43  delete m_slice;
44  }
45 }
46 
47 bool TFCSEnergyAndHitGANV2::is_match_calosample(int calosample) const {
48  if (get_Binning().find(calosample) == get_Binning().cend())
49  return false;
50  if (get_Binning().at(calosample).GetNbinsX() == 1)
51  return false;
52  return true;
53 }
54 
55 unsigned int TFCSEnergyAndHitGANV2::get_nr_of_init(unsigned int bin) const {
56  if (bin >= m_bin_ninit.size())
57  return 0;
58  return m_bin_ninit[bin];
59 }
60 
62  unsigned int ninit) {
63  if (bin >= m_bin_ninit.size()) {
64  m_bin_ninit.resize(bin + 1, 0);
65  m_bin_ninit.shrink_to_fit();
66  }
67  m_bin_ninit[bin] = ninit;
68 }
69 
70 // initialize lwtnn network
72  int const &pid, int const &etaMin,
73  const std::string &FastCaloGANInputFolderName) {
74 
75  // initialize all necessary constants
76  // FIXME eventually all these could be stored in the .json file
77 
79  "Using FastCaloGANInputFolderName: " << FastCaloGANInputFolderName);
80  // get neural net JSON file as an std::istream object
81  const int etaMax = etaMin + 5;
82 
84  set_pdgid(pid);
85  if (pid == 11)
86  add_pdgid(-pid);
87  if (pid == 211)
88  add_pdgid(-pid);
89  set_eta_min(etaMin / 100.0);
90  set_eta_max(etaMax / 100.0);
91  set_eta_nominal((etaMin + etaMax) / 200.0);
92 
93  int pidForXml = pid;
94  if (pid != 22 && pid != 11) {
95  pidForXml = 211;
96  }
97 
98  const int etaMid = (etaMin + etaMax) / 2;
99  m_param.InitialiseFromXML(pidForXml, etaMid, FastCaloGANInputFolderName);
100  m_param.Print();
102  m_slice->Print();
103  return m_slice->LoadGAN();
104 }
105 
107  TFCSSimulationState &simulstate, const TFCSTruthState *,
108  const TFCSExtrapolationState *) const {
109  return std::string(
110  Form("layer=%d", simulstate.getAuxInfo<int>("GANlayer"_FCShash)));
111 }
112 
114  TFCSSimulationState &simulstate, const TFCSTruthState *truth,
115  const TFCSExtrapolationState *extrapol) const {
116  if (!truth) {
117  ATH_MSG_ERROR("Invalid truth pointer");
118  return false;
119  }
120 
121  if (!extrapol) {
122  ATH_MSG_ERROR("Invalid extrapolation pointer");
123  return false;
124  }
125 
126  const int pdgId = truth->pdgid();
127  const float charge = HepPDT::ParticleID(pdgId).charge();
128 
129  float Einit;
130  const float Ekin = truth->Ekin();
131  if (OnlyScaleEnergy())
132  Einit = simulstate.E();
133  else
134  Einit = Ekin;
135 
136  ATH_MSG_VERBOSE("Momentum " << truth->P() << " pdgId " << truth->pdgid());
137  // check that the network exists
138  if (!m_slice->IsGanCorrectlyLoaded()) {
139  ATH_MSG_WARNING("GAN not loaded correctly.");
140  return false;
141  }
142 
144  m_slice->GetNetworkOutputs(truth, extrapol, simulstate);
145  ATH_MSG_VERBOSE("network outputs size: " << outputs.size());
146 
147  const TFCSGANXMLParameters::Binning &binsInLayers = m_param.GetBinning();
148  const auto ganVersion = m_param.GetGANVersion();
149  const TFCSGANEtaSlice::FitResultsPerLayer &fitResults =
150  m_slice->GetFitResults(); // used only if GAN version > 1
151 
152  ATH_MSG_DEBUG("energy voxels size = " << outputs.size());
153 
154  double totalEnergy = 0;
155  for (auto output : outputs) {
156  totalEnergy += output.second;
157  }
158  if (totalEnergy < 0) {
159  ATH_MSG_WARNING("Energy from GAN is negative, skipping particle");
160  return false;
161  }
162 
163  ATH_MSG_VERBOSE("Get binning");
164 
165  simulstate.set_E(0);
166 
167  int vox = 0;
168  for (const auto &[layer, h] : binsInLayers) {
169  // attempt to debug intermittent ci issues described in
170  // https://its.cern.ch/jira/browse/ATLASSIM-7031
171  if (h.IsZombie()) {
172  ATH_MSG_INFO("Got truth state: ");
173  truth->Print();
174 
175  ATH_MSG_INFO("Got extrapolation state: ");
176  extrapol->Print();
177 
178  ATH_MSG_INFO("Got simulation state: ");
179  simulstate.Print();
180 
181  ATH_MSG_INFO("Got GAN XML parameters: ");
182  m_param.Print();
183 
184  ATH_MSG_ERROR("Histogram pointer for layer " << layer << " is broken");
185  return false;
186  }
187 
188  const int xBinNum = h.GetNbinsX();
189  const int yBinNum = h.GetNbinsY();
190  const TAxis *x = h.GetXaxis();
191 
192  // If only one bin in r means layer is empty, no value should be added
193  if (xBinNum == 1) {
194  ATH_MSG_VERBOSE(" Layer "
195  << layer
196  << " has only one bin in r, this means is it not used, "
197  "skipping (this is needed to keep correct "
198  "syncronisation of voxel and layers)");
199  // delete h;
200  continue;
201  }
202 
203  ATH_MSG_VERBOSE(" Getting energy for Layer " << layer);
204 
205  // First fill energies
206  for (int ix = 1; ix <= xBinNum; ++ix) {
207  double binsInAlphaInRBin = GetAlphaBinsForRBin(x, ix, yBinNum);
208  for (int iy = 1; iy <= binsInAlphaInRBin; ++iy) {
209  const double energyInVoxel = outputs.at(std::to_string(vox));
210  ATH_MSG_VERBOSE(" Vox " << vox << " energy " << energyInVoxel
211  << " binx " << ix << " biny " << iy);
212 
213  if (energyInVoxel <= 0) {
214  vox++;
215  continue;
216  }
217 
218  simulstate.add_E(layer, Einit * energyInVoxel);
219  vox++;
220  }
221  }
222  }
223 
224  for (unsigned int ichain = m_bin_start.back(); ichain < size(); ++ichain) {
225  ATH_MSG_DEBUG("now run for all bins: " << chain()[ichain]->GetName());
226  if (simulate_and_retry(chain()[ichain], simulstate, truth, extrapol) !=
227  FCSSuccess) {
228  return FCSFatal;
229  }
230  }
231 
232  vox = 0;
233  for (const auto &[layer, h] : binsInLayers) {
234  const int xBinNum = h.GetNbinsX();
235  const int yBinNum = h.GetNbinsY();
236  const TAxis *x = h.GetXaxis();
237  const TAxis *y = h.GetYaxis();
238 
239  simulstate.setAuxInfo<int>("GANlayer"_FCShash, layer);
241 
242  // If only one bin in r means layer is empty, no value should be added
243  if (xBinNum == 1) {
244  ATH_MSG_VERBOSE(" Layer "
245  << layer
246  << " has only one bin in r, this means is it not used, "
247  "skipping (this is needed to keep correct "
248  "syncronisation of voxel and layers)");
249  // delete h;
250  continue;
251  }
252 
253  if (get_number_of_bins() > 0) {
254  const int bin = get_bin(simulstate, truth, extrapol);
255  if (bin >= 0 && bin < (int)get_number_of_bins()) {
256  for (unsigned int ichain = m_bin_start[bin];
257  ichain < TMath::Min(m_bin_start[bin] + get_nr_of_init(bin),
258  m_bin_start[bin + 1]);
259  ++ichain) {
260  ATH_MSG_DEBUG("for " << get_variable_text(simulstate, truth, extrapol)
261  << " run init " << get_bin_text(bin) << ": "
262  << chain()[ichain]->GetName());
263  if (chain()[ichain]->InheritsFrom(
264  TFCSLateralShapeParametrizationHitBase::Class())) {
266  static_cast<TFCSLateralShapeParametrizationHitBase *>(chain()[ichain]);
267  if (sim->simulate_hit(hit, simulstate, truth, extrapol) !=
268  FCSSuccess) {
269  ATH_MSG_ERROR("error for "
270  << get_variable_text(simulstate, truth, extrapol)
271  << " run init " << get_bin_text(bin) << ": "
272  << chain()[ichain]->GetName());
273  return false;
274  }
275  } else {
276  ATH_MSG_ERROR("for "
277  << get_variable_text(simulstate, truth, extrapol)
278  << " run init " << get_bin_text(bin) << ": "
279  << chain()[ichain]->GetName()
280  << " does not inherit from "
281  "TFCSLateralShapeParametrizationHitBase");
282  return false;
283  }
284  }
285  } else {
286  ATH_MSG_WARNING("nothing to init for "
287  << get_variable_text(simulstate, truth, extrapol)
288  << ": " << get_bin_text(bin));
289  }
290  }
291 
292  int binResolution = 5;
293  if (layer == 1 || layer == 5) {
294  binResolution = 1;
295  }
296 
297  const double center_eta = hit.center_eta();
298  const double center_phi = hit.center_phi();
299  const double center_r = hit.center_r();
300  const double center_z = hit.center_z();
301 
302  ATH_MSG_VERBOSE(" Layer " << layer << " Extrap eta " << center_eta
303  << " phi " << center_phi << " R " << center_r);
304 
305  const float dist000 =
306  TMath::Sqrt(center_r * center_r + center_z * center_z);
307  const float eta_jakobi = TMath::Abs(2.0 * TMath::Exp(-center_eta) /
308  (1.0 + TMath::Exp(-2 * center_eta)));
309 
310  int nHitsAlpha{};
311  int nHitsR{};
312 
313  // Now create hits
314  for (int ix = 1; ix <= xBinNum; ++ix) {
315  const int binsInAlphaInRBin = GetAlphaBinsForRBin(x, ix, yBinNum);
316 
317  // Horrible work around for variable # of bins along alpha direction
318  const int binsToMerge = yBinNum == 32 ? 32 / binsInAlphaInRBin : 1;
319  for (int iy = 1; iy <= binsInAlphaInRBin; ++iy) {
320  const double energyInVoxel = outputs.at(std::to_string(vox));
321  const int lowEdgeIndex = (iy - 1) * binsToMerge + 1;
322 
323  ATH_MSG_VERBOSE(" Vox " << vox << " energy " << energyInVoxel
324  << " binx " << ix << " biny " << iy);
325 
326  if (energyInVoxel <= 0) {
327  vox++;
328  continue;
329  }
330 
331  if (std::abs(pdgId) == 22 || std::abs(pdgId) == 11) {
332  // maximum 10 MeV per hit, equaly distributed in alpha and r
333  int maxHitsInVoxel = energyInVoxel * truth->Ekin() / 10;
334  if (maxHitsInVoxel < 1)
335  maxHitsInVoxel = 1;
336  nHitsAlpha = std::sqrt(maxHitsInVoxel);
337  nHitsR = std::sqrt(maxHitsInVoxel);
338  } else {
339  // One hit per mm along r
340  nHitsR = x->GetBinUpEdge(ix) - x->GetBinLowEdge(ix);
341  if (yBinNum == 1) {
342  // nbins in alpha depend on circumference lenght
343  const double r = x->GetBinUpEdge(ix);
344  nHitsAlpha = ceil(2 * TMath::Pi() * r / binResolution);
345  } else {
346  // d = 2*r*sin (a/2r) this distance at the upper r must be 1mm for
347  // layer 1 or 5, 5mm otherwise.
348  const double angle = y->GetBinUpEdge(iy) - y->GetBinLowEdge(iy);
349  const double r = x->GetBinUpEdge(ix);
350  const double d = 2 * r * sin(angle / 2 * r);
351  nHitsAlpha = ceil(d / binResolution);
352  }
353 
354  if (layer != 1 && layer != 5) {
355  // For layers that are not EMB1 or EMEC1 use a maximum of 10 hits
356  // per direction, a higher granularity is needed for the other
357  // layers
358  const int maxNhits = 10;
359  nHitsAlpha = std::min(maxNhits, std::max(1, nHitsAlpha));
360  nHitsR = std::min(maxNhits, std::max(1, nHitsR));
361  }
362  }
363 
364  for (int ir = 0; ir < nHitsR; ++ir) {
365  double r =
366  x->GetBinLowEdge(ix) + x->GetBinWidth(ix) / (nHitsR + 1) * ir;
367 
368  for (int ialpha = 1; ialpha <= nHitsAlpha; ++ialpha) {
369  if (ganVersion > 1) {
370  if (fitResults.at(layer)[ix - 1] != 0) {
371  int tries = 0;
372  double a = CLHEP::RandFlat::shoot(simulstate.randomEngine(),
373  x->GetBinLowEdge(ix),
374  x->GetBinUpEdge(ix));
375  double rand_r =
376  log((a - x->GetBinLowEdge(ix)) / (x->GetBinWidth(ix))) /
377  fitResults.at(layer)[ix - 1];
378  while ((rand_r < x->GetBinLowEdge(ix) ||
379  rand_r > x->GetBinUpEdge(ix)) &&
380  tries < 100) {
381  a = CLHEP::RandFlat::shoot(simulstate.randomEngine(),
382  x->GetBinLowEdge(ix),
383  x->GetBinUpEdge(ix));
384  rand_r =
385  log((a - x->GetBinLowEdge(ix)) / (x->GetBinWidth(ix))) /
386  fitResults.at(layer)[ix - 1];
387  tries++;
388  }
389  if (tries >= 100) {
390  ATH_MSG_VERBOSE(" Too many tries for bin ["
391  << x->GetBinLowEdge(ix) << "-"
392  << x->GetBinUpEdge(ix) << "] having slope "
393  << fitResults.at(layer)[ix - 1]
394  << " will use grid (old method)");
395  } else {
396  r = rand_r;
397  }
398  }
399  }
400 
401  double alpha;
402  if (binsInAlphaInRBin == 1) {
403  alpha = CLHEP::RandFlat::shoot(simulstate.randomEngine(),
404  -TMath::Pi(), TMath::Pi());
405  } else {
406  alpha =
407  y->GetBinLowEdge(lowEdgeIndex) +
408  y->GetBinWidth(iy) * binsToMerge / (nHitsAlpha + 1) * ialpha;
409 
410  if (m_param.IsSymmetrisedAlpha()) {
411  if (CLHEP::RandFlat::shoot(simulstate.randomEngine(),
412  -TMath::Pi(), TMath::Pi()) < 0) {
413  alpha = -alpha;
414  }
415  }
416  }
417 
418  hit.reset();
419  hit.E() = Einit * energyInVoxel / (nHitsAlpha * nHitsR);
420 
421  if (layer <= 20) {
422  float delta_eta_mm = r * cos(alpha);
423  float delta_phi_mm = r * sin(alpha);
424 
425  ATH_MSG_VERBOSE("delta_eta_mm " << delta_eta_mm
426  << " delta_phi_mm "
427  << delta_phi_mm);
428 
429  // Particles with negative eta are expected to have the same shape
430  // as those with positive eta after transformation: delta_eta -->
431  // -delta_eta
432  if (center_eta < 0.)
433  delta_eta_mm = -delta_eta_mm;
434  // We derive the shower shapes for electrons and positively
435  // charged hadrons. Particle with the opposite charge are expected
436  // to have the same shower shape after the transformation:
437  // delta_phi --> -delta_phi
438  if ((charge < 0. && pdgId != 11) || pdgId == -11)
439  delta_phi_mm = -delta_phi_mm;
440 
441  const float delta_eta = delta_eta_mm / eta_jakobi / dist000;
442  const float delta_phi = delta_phi_mm / center_r;
443 
444  hit.eta() = center_eta + delta_eta;
445  hit.phi() = TVector2::Phi_mpi_pi(center_phi + delta_phi);
446 
447  ATH_MSG_VERBOSE(" Hit eta " << hit.eta() << " phi " << hit.phi()
448  << " layer " << layer);
449  } else { // FCAL is in (x,y,z)
450  const float hit_r = r * cos(alpha) + center_r;
451  float delta_phi = r * sin(alpha) / center_r;
452  // We derive the shower shapes for electrons and positively
453  // charged hadrons. Particle with the opposite charge are expected
454  // to have the same shower shape after the transformation:
455  // delta_phi --> -delta_phi
456  if ((charge < 0. && pdgId != 11) || pdgId == -11)
457  delta_phi = -delta_phi;
458  const float hit_phi =
459  TVector2::Phi_mpi_pi(center_phi + delta_phi);
460  hit.x() = hit_r * cos(hit_phi);
461  hit.y() = hit_r * sin(hit_phi);
462  hit.z() = center_z;
463  ATH_MSG_VERBOSE(" Hit x " << hit.x() << " y " << hit.y()
464  << " layer " << layer);
465  }
466 
467  if (get_number_of_bins() > 0) {
468  const int bin = get_bin(simulstate, truth, extrapol);
469  if (bin >= 0 && bin < (int)get_number_of_bins()) {
470  for (unsigned int ichain =
472  ichain < m_bin_start[bin + 1]; ++ichain) {
474  "for " << get_variable_text(simulstate, truth, extrapol)
475  << " run " << get_bin_text(bin) << ": "
476  << chain()[ichain]->GetName());
477  if (chain()[ichain]->InheritsFrom(
478  TFCSLateralShapeParametrizationHitBase::Class())) {
481  *>(chain()[ichain]);
482  if (sim->simulate_hit(hit, simulstate, truth, extrapol) !=
483  FCSSuccess) {
485  "error for "
486  << get_variable_text(simulstate, truth, extrapol)
487  << " run init " << get_bin_text(bin) << ": "
488  << chain()[ichain]->GetName());
489  return false;
490  }
491  } else {
493  "for " << get_variable_text(simulstate, truth, extrapol)
494  << " run init " << get_bin_text(bin) << ": "
495  << chain()[ichain]->GetName()
496  << " does not inherit from "
497  "TFCSLateralShapeParametrizationHitBase");
498  return false;
499  }
500  }
501  } else {
503  "nothing to do for "
504  << get_variable_text(simulstate, truth, extrapol) << ": "
505  << get_bin_text(bin));
506  }
507  } else {
508  ATH_MSG_WARNING("no bins defined, is this intended?");
509  }
510  }
511  }
512  vox++;
513  }
514  }
515 
516  ATH_MSG_VERBOSE("Number of voxels " << vox);
517  ATH_MSG_VERBOSE("Done layer " << layer);
518  }
519 
520  if (simulstate.E() > std::numeric_limits<double>::epsilon()) {
521  for (int ilayer = 0; ilayer < CaloCell_ID_FCS::MaxSample; ++ilayer) {
522  simulstate.set_Efrac(ilayer, simulstate.E(ilayer) / simulstate.E());
523  }
524  }
525 
526  ATH_MSG_VERBOSE("Done particle");
527  return true;
528 }
529 
531  TFCSSimulationState &simulstate, const TFCSTruthState *truth,
532  const TFCSExtrapolationState *extrapol) const {
533  for (unsigned int ichain = 0; ichain < m_bin_start[0]; ++ichain) {
534  ATH_MSG_DEBUG("now run for all bins: " << chain()[ichain]->GetName());
535  if (simulate_and_retry(chain()[ichain], simulstate, truth, extrapol) !=
536  FCSSuccess) {
537  return FCSFatal;
538  }
539  }
540 
541  ATH_MSG_VERBOSE("Fill Energies");
542  if (!fillEnergy(simulstate, truth, extrapol)) {
543  ATH_MSG_WARNING("Could not fill energies ");
544  // bail out but do not stop the job
545  return FCSSuccess;
546  }
547 
548  return FCSSuccess;
549 }
550 
551 void TFCSEnergyAndHitGANV2::Print(Option_t *option) const {
553  TString opt(option);
554  const bool shortprint = opt.Index("short") >= 0;
555  const bool longprint =
556  msgLvl(MSG::DEBUG) || (msgLvl(MSG::INFO) && !shortprint);
557  TString optprint = opt;
558  optprint.ReplaceAll("short", "");
559 
560  TString prefix = "- ";
561  for (unsigned int ichain = 0; ichain < size(); ++ichain) {
562  if (ichain == 0 && ichain != m_bin_start.front()) {
563  prefix = "> ";
564  if (longprint)
565  ATH_MSG_INFO(optprint << prefix << "Run for all bins");
566  }
567  for (unsigned int ibin = 0; ibin < get_number_of_bins(); ++ibin) {
568  if (ichain == m_bin_start[ibin]) {
569  if (ibin < get_number_of_bins() - 1)
570  if (ichain == m_bin_start[ibin + 1])
571  continue;
572  prefix = Form("%-2d", ibin);
573  if (longprint)
574  ATH_MSG_INFO(optprint << prefix << "Run for " << get_bin_text(ibin));
575  }
576  }
577  if (ichain == m_bin_start.back()) {
578  prefix = "< ";
579  if (longprint)
580  ATH_MSG_INFO(optprint << prefix << "Run for all bins");
581  }
582  chain()[ichain]->Print(opt + prefix);
583  }
584 }
585 
587  const TFCSTruthState *truth,
590  ATH_MSG_NOCLASS(logger, "Start lwtnn test" << std::endl);
591  std::string path =
592  "/eos/atlas/atlascerngroupdisk/proj-simul/AF3_Run3/"
593  "InputsToBigParamFiles/FastCaloGANWeightsVer02/";
594  test_path(path, simulstate, truth, extrapol, "lwtnn");
595 
596  ATH_MSG_NOCLASS(logger, "Start onnx test" << std::endl);
597  path =
598  "/eos/atlas/atlascerngroupdisk/proj-simul/AF3_Run3/"
599  "InputsToBigParamFiles/FastCaloGANWeightsONNXVer08/";
600  test_path(path, simulstate, truth, extrapol, "onnx");
601  ATH_MSG_NOCLASS(logger, "Finish all tests" << std::endl);
602 }
603 
604 void TFCSEnergyAndHitGANV2::test_path(const std::string &path,
605  TFCSSimulationState *simulstate,
606  const TFCSTruthState *truth,
608  const std::string &outputname, int pid) {
610  ATH_MSG_NOCLASS(logger, "Running test on " << path << std::endl);
611  if (!simulstate) {
612  simulstate = new TFCSSimulationState();
613 #if defined(__FastCaloSimStandAlone__)
614  simulstate->setRandomEngine(new CLHEP::TRandomEngine());
615 #else
616  simulstate->setRandomEngine(new CLHEP::RanluxEngine());
617 #endif
618  }
619  if (!truth) {
620  ATH_MSG_NOCLASS(logger, "New particle");
622  t->SetPtEtaPhiM(65536, 0, 0, ParticleConstants::chargedPionMassInMeV);
623  t->set_pdgid(pid);
624  truth = t;
625  }
626  if (!extrapol) {
628  e->set_IDCaloBoundary_eta(truth->Eta());
629  for (int i = 0; i < 24; ++i) {
630  e->set_eta(i, TFCSExtrapolationState::SUBPOS_ENT, truth->Eta());
631  e->set_eta(i, TFCSExtrapolationState::SUBPOS_EXT, truth->Eta());
632  e->set_eta(i, TFCSExtrapolationState::SUBPOS_MID, truth->Eta());
633  e->set_phi(i, TFCSExtrapolationState::SUBPOS_ENT, 0);
634  e->set_phi(i, TFCSExtrapolationState::SUBPOS_EXT, 0);
635  e->set_phi(i, TFCSExtrapolationState::SUBPOS_MID, 0);
636  e->set_r(i, TFCSExtrapolationState::SUBPOS_ENT, 1500 + i * 10);
637  e->set_r(i, TFCSExtrapolationState::SUBPOS_EXT, 1510 + i * 10);
638  e->set_r(i, TFCSExtrapolationState::SUBPOS_MID, 1505 + i * 10);
639  e->set_z(i, TFCSExtrapolationState::SUBPOS_ENT, 3500 + i * 10);
640  e->set_z(i, TFCSExtrapolationState::SUBPOS_EXT, 3510 + i * 10);
641  e->set_z(i, TFCSExtrapolationState::SUBPOS_MID, 3505 + i * 10);
642  }
643  extrapol = e;
644  }
645 
646  TFCSEnergyAndHitGANV2 GAN("GAN", "GAN");
647  GAN.setLevel(MSG::INFO);
648  const int etaMin = 20;
649  const int etaMax = etaMin + 5;
650  ATH_MSG_NOCLASS(logger, "Initialize Networks");
652  for (int i = 0; i < 24; ++i)
653  if (GAN.is_match_calosample(i)) {
655  Form("center%d", i), Form("center layer %d", i));
656  c->set_calosample(i);
657  c->setExtrapWeight(0.5);
658  c->setLevel(MSG::INFO);
659  c->set_pdgid(pid);
660  if (pid == 11)
661  c->add_pdgid(-pid);
662  if (pid == 211)
663  c->add_pdgid(-pid);
664  c->set_eta_min(etaMin / 100.0);
665  c->set_eta_max(etaMax / 100.0);
666  c->set_eta_nominal((etaMin + etaMax) / 200.0);
667 
668  GAN.push_back_in_bin(c, i);
669  GAN.set_nr_of_init(i, 1);
670  }
671 
672  GAN.Print();
673 
674  ATH_MSG_NOCLASS(logger, "Writing GAN to " << outputname);
675  const std::string outname = "FCSGANtest_" + outputname + ".root";
676  TFile *fGAN = TFile::Open(outname.c_str(), "recreate");
677  fGAN->cd();
678  // GAN.Write();
679  fGAN->WriteObjectAny(&GAN, "TFCSEnergyAndHitGANV2", "GAN");
680 
681  fGAN->ls();
682  fGAN->Close();
683 
684  ATH_MSG_NOCLASS(logger, "Open " << outname);
685  fGAN = TFile::Open(outname.c_str());
686  TFCSEnergyAndHitGANV2 *GAN2 = (TFCSEnergyAndHitGANV2 *)(fGAN->Get("GAN"));
687  GAN2->setLevel(MSG::INFO);
688  GAN2->Print();
689 
690  ATH_MSG_NOCLASS(logger, "Before running GAN2->simulate()");
691  GAN2->simulate(*simulstate, truth, extrapol);
692  simulstate->Print();
693 }
694 
696  if (bins < 4)
697  return 4;
698  else if (bins < 8)
699  return 8;
700  else if (bins < 16)
701  return 16;
702  else
703  return 32;
704 }
705 
707  int yBinNum) const {
708  double binsInAlphaInRBin = yBinNum;
709  if (yBinNum == 32) {
710  ATH_MSG_DEBUG("yBinNum is special value 32");
711  const double widthX = x->GetBinWidth(ix);
712  const double radious = x->GetBinCenter(ix);
713  double circumference = radious * 2 * TMath::Pi();
714  if (m_param.IsSymmetrisedAlpha()) {
715  circumference = radious * TMath::Pi();
716  }
717 
718  const double bins = circumference / widthX;
719  binsInAlphaInRBin = GetBinsInFours(bins);
720  ATH_MSG_DEBUG("Bin in alpha: " << binsInAlphaInRBin << " for r bin: " << ix
721  << " (" << x->GetBinLowEdge(ix) << "-"
722  << x->GetBinUpEdge(ix) << ")");
723  }
724  return binsInAlphaInRBin;
725 }
TFCSCenterPositionCalculation.h
TFCSGANXMLParameters::Print
void Print() const
Definition: TFCSGANXMLParameters.cxx:124
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
beamspotman.r
def r
Definition: beamspotman.py:674
FCSReturnCode
FCSReturnCode
Base class for all FastCaloSim parametrizations Functionality in derivde classes is provided through ...
Definition: TFCSParametrizationBase.h:41
TFCSGANEtaSlice::LoadGAN
bool LoadGAN()
Definition: TFCSGANEtaSlice.cxx:71
TFCSLateralShapeParametrizationHitBase::simulate_hit
virtual FCSReturnCode simulate_hit(Hit &hit, TFCSSimulationState &simulstate, const TFCSTruthState *truth, const TFCSExtrapolationState *extrapol)
simulated one hit position with some energy.
Definition: TFCSLateralShapeParametrizationHitBase.cxx:50
ISF_FCS::MLogging
Cut down AthMessaging.
Definition: MLogging.h:176
TFCSSimulationState::getAuxInfo
const T getAuxInfo(std::uint32_t index) const
Definition: TFCSSimulationState.h:161
TFCSParametrizationBinnedChain::get_bin_text
virtual const std::string get_bin_text(int bin) const
print the range of a bin; for bin -1, print the allowed range
Definition: TFCSParametrizationBinnedChain.cxx:60
TFCSEnergyAndHitGANV2::set_GANfreemem
void set_GANfreemem()
Definition: TFCSEnergyAndHitGANV2.h:40
TFCSEnergyAndHitGANV2
Definition: TFCSEnergyAndHitGANV2.h:20
add-xsec-uncert-quadrature-N.alpha
alpha
Definition: add-xsec-uncert-quadrature-N.py:110
athena.path
path
python interpreter configuration --------------------------------------—
Definition: athena.py:128
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
TFCSEnergyAndHitGANV2::get_Binning
const TFCSGANXMLParameters::Binning & get_Binning() const
Definition: TFCSEnergyAndHitGANV2.h:66
make_coralServer_rep.opt
opt
Definition: make_coralServer_rep.py:19
TFCSEnergyAndHitGANV2::initializeNetwork
bool initializeNetwork(int const &pid, int const &etaMin, const std::string &FastCaloGANInputFolderName)
Definition: TFCSEnergyAndHitGANV2.cxx:71
TFCSGANXMLParameters::InitialiseFromXML
void InitialiseFromXML(int pid, int etaMid, const std::string &FastCaloGANInputFolderName)
Definition: TFCSGANXMLParameters.cxx:16
hist_file_dump.d
d
Definition: hist_file_dump.py:142
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
TFCSGANEtaSlice::NetworkOutputs
std::map< std::string, double > NetworkOutputs
Definition: TFCSGANEtaSlice.h:40
TFCSParametrization::add_pdgid
virtual void add_pdgid(int id)
Definition: TFCSParametrization.cxx:35
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
TFCSEnergyAndHitGANV2::GetBinsInFours
static int GetBinsInFours(double const &bins)
Definition: TFCSEnergyAndHitGANV2.cxx:695
TFCSCenterPositionCalculation
Definition: TFCSCenterPositionCalculation.h:11
TFCSParametrizationBinnedChain::push_back_in_bin
virtual void push_back_in_bin(TFCSParametrizationBase *param, unsigned int bin)
Definition: TFCSParametrizationBinnedChain.cxx:31
bin
Definition: BinsDiffFromStripMedian.h:43
xAOD::etaMax
etaMax
Definition: HIEventShape_v2.cxx:46
TFCSGANXMLParameters::Binning
std::map< int, TH2D > Binning
Definition: TFCSGANXMLParameters.h:23
TFCSExtrapolationState::SUBPOS_ENT
@ SUBPOS_ENT
Definition: TFCSExtrapolationState.h:21
TFCSEnergyAndHitGANV2::OnlyScaleEnergy
bool OnlyScaleEnergy() const
Definition: TFCSEnergyAndHitGANV2.h:50
TFCSExtrapolationState
Definition: TFCSExtrapolationState.h:13
TFCSSimulationState::randomEngine
CLHEP::HepRandomEngine * randomEngine()
Definition: TFCSSimulationState.h:36
TFCSLateralShapeParametrizationHitBase::Hit
Definition: TFCSLateralShapeParametrizationHitBase.h:42
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
TFCSExtrapolationState::SUBPOS_MID
@ SUBPOS_MID
Definition: TFCSExtrapolationState.h:20
RunActsMaterialValidation.extrapol
extrapol
Definition: RunActsMaterialValidation.py:91
TFCSEnergyAndHitGANV2::set_nr_of_init
void set_nr_of_init(unsigned int bin, unsigned int ninit)
Definition: TFCSEnergyAndHitGANV2.cxx:61
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
TFCSEnergyAndHitGANV2.h
TFCSTruthState::Print
void Print(Option_t *option="") const
Definition: TFCSTruthState.cxx:21
Phi_mpi_pi
__HOSTDEV__ double Phi_mpi_pi(double)
Definition: GeoRegion.cxx:7
x
#define x
TFCSExtrapolationState::SUBPOS_EXT
@ SUBPOS_EXT
Definition: TFCSExtrapolationState.h:22
TFCSParametrization::set_eta_max
virtual void set_eta_max(double max)
Definition: TFCSParametrization.cxx:53
TFCSLateralShapeParametrizationHitBase
Definition: TFCSLateralShapeParametrizationHitBase.h:13
TFCSParametrization::set_pdgid
virtual void set_pdgid(int id)
Definition: TFCSParametrization.cxx:28
ParticleConstants::PDG2011::chargedPionMassInMeV
constexpr double chargedPionMassInMeV
the mass of the charged pion (in MeV)
Definition: ParticleConstants.h:41
TFCSParametrizationBinnedChain::get_number_of_bins
virtual unsigned int get_number_of_bins() const
Definition: TFCSParametrizationBinnedChain.h:22
TFCSParametrizationBinnedChain
Definition: TFCSParametrizationBinnedChain.h:10
TFCSSimulationState::add_E
void add_E(int sample, double Esample)
Definition: TFCSSimulationState.h:53
TFCSLateralShapeParametrizationHitBase.h
CaloCell_ID_FCS::MaxSample
@ MaxSample
Definition: FastCaloSim_CaloCell_ID.h:47
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
ParticleConstants.h
lumiFormat.i
int i
Definition: lumiFormat.py:85
TFCSParametrizationChain::chain
const Chain_t & chain() const
Definition: TFCSParametrizationChain.h:53
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
angle
double angle(const GeoTrf::Vector2D &a, const GeoTrf::Vector2D &b)
Definition: TRTDetectorFactory_Full.cxx:68
TFCSTruthState::Ekin
double Ekin() const
Definition: TFCSTruthState.h:26
TRT::Hit::layer
@ layer
Definition: HitInfo.h:79
covarianceTool.title
title
Definition: covarianceTool.py:542
TFCSEnergyAndHitGANV2::m_slice
TFCSGANEtaSlice * m_slice
Definition: TFCSEnergyAndHitGANV2.h:109
ParticleGun_EoverP_Config.pid
pid
Definition: ParticleGun_EoverP_Config.py:62
checkCorrelInHIST.prefix
dictionary prefix
Definition: checkCorrelInHIST.py:391
TFCSGANEtaSlice::Print
void Print() const
Definition: TFCSGANEtaSlice.cxx:258
jobOptions.ParticleID
ParticleID
Definition: jobOptions.decayer.py:85
TFCSEnergyAndHitGANV2::~TFCSEnergyAndHitGANV2
virtual ~TFCSEnergyAndHitGANV2()
Definition: TFCSEnergyAndHitGANV2.cxx:41
TFCSGANEtaSlice
Definition: TFCSGANEtaSlice.h:30
TFCSEnergyAndHitGANV2::GetAlphaBinsForRBin
int GetAlphaBinsForRBin(const TAxis *x, int ix, int yBinNum) const
Definition: TFCSEnergyAndHitGANV2.cxx:706
TFCSParametrizationBase::Print
void Print(Option_t *option="") const
Print object information.
Definition: TFCSParametrizationBase.cxx:52
TFCSGANEtaSlice::FitResultsPerLayer
std::map< int, std::vector< double > > FitResultsPerLayer
Definition: TFCSGANEtaSlice.h:37
plotting.yearwise_luminosity_vs_mu.bins
bins
Definition: yearwise_luminosity_vs_mu.py:30
TFCSSimulationState::Print
void Print(Option_t *option="") const
Definition: TFCSSimulationState.cxx:40
FCSFatal
@ FCSFatal
Definition: TFCSParametrizationBase.h:41
TFCSGANEtaSlice::IsGanCorrectlyLoaded
bool IsGanCorrectlyLoaded() const
Definition: TFCSGANEtaSlice.cxx:58
FCSSuccess
@ FCSSuccess
Definition: TFCSParametrizationBase.h:41
TFCSGANEtaSlice::GetFitResults
const FitResultsPerLayer & GetFitResults() const
Definition: TFCSGANEtaSlice.h:51
python.CreateTierZeroArgdict.outputs
outputs
Definition: CreateTierZeroArgdict.py:189
merge.output
output
Definition: merge.py:16
beamspotman.outname
outname
Definition: beamspotman.py:412
TFCSParametrizationBase::reset_match_all_pdgid
virtual void reset_match_all_pdgid()
Definition: TFCSParametrizationBase.h:84
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
TFCSParametrization::set_eta_nominal
virtual void set_eta_nominal(double min)
Definition: TFCSParametrization.cxx:47
ActsTrk::to_string
std::string to_string(const DetectorType &type)
Definition: GeometryDefs.h:34
ISF_FCS::MLogging::setLevel
virtual void setLevel(MSG::Level lvl)
Update outputlevel.
Definition: MLogging.cxx:105
TFCSSimulationState::setAuxInfo
void setAuxInfo(std::uint32_t index, const T &val)
Definition: TFCSSimulationState.h:168
plotBeamSpotVxVal.bin
int bin
Definition: plotBeamSpotVxVal.py:82
TFCSSimulationState::set_E
void set_E(int sample, double Esample)
Definition: TFCSSimulationState.h:48
charge
double charge(const T &p)
Definition: AtlasPID.h:986
runIDPVM.pdgId
pdgId
Definition: runIDPVM.py:91
TFCSTruthState::pdgid
int pdgid() const
Definition: TFCSTruthState.h:25
TFCSEnergyAndHitGANV2::unit_test
static void unit_test(TFCSSimulationState *simulstate=nullptr, const TFCSTruthState *truth=nullptr, const TFCSExtrapolationState *extrapol=nullptr)
Definition: TFCSEnergyAndHitGANV2.cxx:586
TFCSGANXMLParameters::GetGANVersion
int GetGANVersion() const
Definition: TFCSGANXMLParameters.h:35
TFCSEnergyAndHitGANV2::simulate
virtual FCSReturnCode simulate(TFCSSimulationState &simulstate, const TFCSTruthState *truth, const TFCSExtrapolationState *extrapol) const override
Method in all derived classes to do some simulation.
Definition: TFCSEnergyAndHitGANV2.cxx:530
TFCSGANEtaSlice::GetNetworkOutputs
NetworkOutputs GetNetworkOutputs(const TFCSTruthState *truth, const TFCSExtrapolationState *extrapol, TFCSSimulationState simulstate) const
Definition: TFCSGANEtaSlice.cxx:174
TFCSSimulationState::setRandomEngine
void setRandomEngine(CLHEP::HepRandomEngine *engine)
Definition: TFCSSimulationState.h:37
TFCSEnergyAndHitGANV2::m_bin_ninit
std::vector< int > m_bin_ninit
Definition: TFCSEnergyAndHitGANV2.h:102
ir
int ir
counter of the current depth
Definition: fastadd.cxx:49
a
TList * a
Definition: liststreamerinfos.cxx:10
LArCellBinning.etaMin
etaMin
Definition: LArCellBinning.py:84
TFCSTruthState.h
TFCSEnergyAndHitGANV2::m_param
TFCSGANXMLParameters m_param
Definition: TFCSEnergyAndHitGANV2.h:110
ATH_MSG_NOCLASS
#define ATH_MSG_NOCLASS(logger_name, x)
Definition: MLogging.h:52
y
#define y
h
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
python.Constants.INFO
int INFO
Definition: Control/AthenaCommon/python/Constants.py:15
TFCSEnergyAndHitGANV2::is_match_calosample
virtual bool is_match_calosample(int calosample) const override
Definition: TFCSEnergyAndHitGANV2.cxx:47
TFCSExtrapolationState.h
TFCSEnergyAndHitGANV2::get_bin
virtual int get_bin(TFCSSimulationState &simulstate, const TFCSTruthState *, const TFCSExtrapolationState *) const override
use the layer to be done as binning of the GAN chain
Definition: TFCSEnergyAndHitGANV2.h:55
eFEXNTuple.delta_phi
def delta_phi(phi1, phi2)
Definition: eFEXNTuple.py:15
DEBUG
#define DEBUG
Definition: page_access.h:11
TFCSParametrizationChain::simulate_and_retry
FCSReturnCode simulate_and_retry(TFCSParametrizationBase *parametrization, TFCSSimulationState &simulstate, const TFCSTruthState *truth, const TFCSExtrapolationState *extrapol) const
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
TFCSParametrization::set_eta_min
virtual void set_eta_min(double min)
Definition: TFCSParametrization.cxx:51
TFCSSimulationState::E
double E() const
Definition: TFCSSimulationState.h:42
TFCSSimulationState.h
TFCSEnergyAndHitGANV2::fillEnergy
bool fillEnergy(TFCSSimulationState &simulstate, const TFCSTruthState *truth, const TFCSExtrapolationState *extrapol) const
Definition: TFCSEnergyAndHitGANV2.cxx:113
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
TFCSEnergyAndHitGANV2::get_nr_of_init
unsigned int get_nr_of_init(unsigned int bin) const
Definition: TFCSEnergyAndHitGANV2.cxx:55
TFCSTruthState
Definition: TFCSTruthState.h:13
TFCSEnergyAndHitGANV2::TFCSEnergyAndHitGANV2
TFCSEnergyAndHitGANV2(const char *name=nullptr, const char *title=nullptr)
Definition: TFCSEnergyAndHitGANV2.cxx:35
python.compressB64.c
def c
Definition: compressB64.py:93
TFCSEnergyAndHitGANV2::test_path
static void test_path(const std::string &path, TFCSSimulationState *simulstate=nullptr, const TFCSTruthState *truth=nullptr, const TFCSExtrapolationState *extrapol=nullptr, const std::string &outputname="unnamed", int pid=211)
Definition: TFCSEnergyAndHitGANV2.cxx:604
TFCSEnergyAndHitGANV2::Print
virtual void Print(Option_t *option="") const override
Definition: TFCSEnergyAndHitGANV2.cxx:551
TFCSSimulationState
Definition: TFCSSimulationState.h:32
TFCSGANXMLParameters::IsSymmetrisedAlpha
bool IsSymmetrisedAlpha() const
Definition: TFCSGANXMLParameters.h:36
python.iconfTool.gui.pad.logger
logger
Definition: pad.py:14
ISF_FCS::MLogging::msgLvl
bool msgLvl(const MSG::Level lvl) const
Check whether the logging system is active at the provided verbosity level.
Definition: MLogging.h:222
TFCSParametrizationBinnedChain::m_bin_start
std::vector< unsigned int > m_bin_start
Contains the index where the TFCSParametrizationBase* instances to run for a given bin start.
Definition: TFCSParametrizationBinnedChain.h:53
TFCSGANXMLParameters::GetBinning
const Binning & GetBinning() const
Definition: TFCSGANXMLParameters.h:33
TFCSSimulationState::set_Efrac
void set_Efrac(int sample, double Efracsample)
Definition: TFCSSimulationState.h:49
TFCSEnergyAndHitGANV2::get_variable_text
virtual const std::string get_variable_text(TFCSSimulationState &simulstate, const TFCSTruthState *, const TFCSExtrapolationState *) const override
Definition: TFCSEnergyAndHitGANV2.cxx:106
TFCSParametrizationChain::size
virtual unsigned int size() const override
Some derived classes have daughter instances of TFCSParametrizationBase objects The size() and operat...
Definition: TFCSParametrizationChain.h:41