34#include "CaloDetDescr/CaloDetDescrElement.h"
43#include "lwtnn/parse_json.hh"
47#include "CLHEP/Random/RandGaussZiggurat.h"
48#include "CLHEP/Random/RandFlat.h"
49#include "CLHEP/Units/SystemOfUnits.h"
64 declareProperty(
"ParamsInputFilename" ,
m_paramsFilename=
"DNNCaloSim/DNNCaloSim_GAN_nn_v0.json",
" lwtnn json output describing the trained network");
65 declareProperty(
"ParamsInputArchitecture" ,
m_paramsInputArchitecture=
"GANv0",
"tag describing additional parameters necessary for the network");
88 ATH_MSG_ERROR(
"Could not get random number engine from RandomNumberService. Abort.");
89 return StatusCode::FAILURE;
100 return StatusCode::FAILURE;
109 m_caloGeo = std::make_unique<CaloGeometryFromCaloDDM>();
111 if(!
m_caloGeo->LoadFCalChannelMapFromFCalDDM(fcalManager) )
ATH_MSG_FATAL(
"Found inconsistency between FCal_Channel map and GEO file. Please, check if they are configured properly.");
118 return StatusCode::FAILURE;
127 return StatusCode::FAILURE;
134 return StatusCode::SUCCESS;
145 if (inputFile.empty()){
147 return StatusCode::FAILURE;
153 std::ifstream input(inputFile);
155 m_graph=std::make_unique<lwt::LightweightGraph>(lwt::parse_json_graph(input) );
158 return StatusCode::FAILURE;
180 return StatusCode::FAILURE;
183 return StatusCode::SUCCESS;
190 return StatusCode::SUCCESS;
195 const EventContext& ctx = Gaudi::Hive::currentContext();
204 return StatusCode::FAILURE;
212 for (; itrTool != endTool; ++itrTool)
214 std::string chronoName=this->name()+
"_"+ itrTool->name();
225 return StatusCode::FAILURE;
237 const int ntrial=100;
238 ATH_MSG_INFO (
"Trial window building on " << ntrial <<
" dummy eta phi " );
239 for (
int i=0 ; i< ntrial ; i++){
240 const double eta = CLHEP::RandFlat::shoot(testsimulstate.
randomEngine(), 0.2, 0.25);
241 const double phi = CLHEP::RandFlat::shoot(testsimulstate.
randomEngine(), -CLHEP::pi , CLHEP::pi);
248 ATH_MSG_INFO (
"End of trial window building on " << ntrial <<
" dummy eta phi " );
255 return StatusCode::SUCCESS;
260 const EventContext& ctx = Gaudi::Hive::currentContext();
268 for (; itrTool != endTool; ++itrTool)
280 return StatusCode::SUCCESS;
288 float aEtaRaw =
a->caloDDE()->eta_raw();
289 float bEtaRaw = b->caloDDE()->eta_raw();
291 float aPhiRaw =
a->caloDDE()->phi_raw();
292 float bPhiRaw = b->caloDDE()->phi_raw();
294 if ((aSampling) < (bSampling))
296 else if ((aSampling) > (bSampling))
299 if ((aEtaRaw) < (bEtaRaw))
301 else if ((aEtaRaw) > (bEtaRaw))
312 float aEtaRaw =
a->caloDDE()->eta_raw();
313 float bEtaRaw = b->caloDDE()->eta_raw();
315 float aPhiRaw =
a->caloDDE()->phi_raw();
316 float bPhiRaw = b->caloDDE()->phi_raw();
318 if ((aSampling) < (bSampling))
320 else if ((aSampling) > (bSampling))
323 if ((aEtaRaw) < (bEtaRaw))
325 else if ((aEtaRaw) > (bEtaRaw))
335 ATH_MSG_VERBOSE(
"NEW PARTICLE! DNNCaloSimSvc called with ISFParticle: " << isfp);
339 if(isfp.
ekin() < 10) {
340 ATH_MSG_VERBOSE(
"Skipping particle with Ekin: " << isfp.
ekin() <<
" MeV. Below the 10 MeV threshold.");
341 return StatusCode::SUCCESS;
351 return StatusCode::SUCCESS;
364 windowCell->addEnergy(trueEnergy * outputs[std::to_string(itr)]);
368 <<
" phi_raw " << windowCell->caloDDE()->phi_raw()
369 <<
" sampling " << windowCell->caloDDE()->getSampling()
370 <<
" energy " << windowCell->energy());
377 return StatusCode::SUCCESS;
392 trueEnergy = isfp.
ekin();
407 for (
int isubpos=0; isubpos< 3 ; isubpos++){
410 " isubpos=" << isubpos <<
411 " OK=" << extrapol.OK(isam,isubpos) <<
412 " eta=" << extrapol.eta(isam,isubpos) <<
413 " phi=" << extrapol.phi(isam,isubpos) <<
414 " r=" << extrapol.r(isam,isubpos) );
423 double etaExtrap=-999.;
424 double phiExtrap=-999.;
425 if (extrapol.eta(isam,isubpos)) {
426 etaExtrap=extrapol.eta(isam,isubpos);
427 phiExtrap=extrapol.phi(isam,isubpos);
431 " isubpos=" << isubpos <<
432 " eta=" << etaExtrap <<
433 " phi=" << phiExtrap );
441 return StatusCode::FAILURE;
450 double randGaussz = 0.;
455 double etaRawImpactCell=impactCellDDE->
eta_raw();
456 double phiRawImpactCell=impactCellDDE->
phi_raw();
459 <<
" phi_index " << impactPhiIndex
462 int pconf = impactPhiIndex % 4 ;
463 int econf = (impactEtaIndex + 1) % 2 ;
465 double riImpactEta = (etaExtrap - etaRawImpactCell);
482 randGaussz = CLHEP::RandGaussZiggurat::shoot(simulstate.
randomEngine(), 0., 1.);
483 inputs[
"Z"].insert ( std::pair<std::string,double>(std::to_string(i), randGaussz) );
490 for (
int i = 0; i< 4; i ++)
493 inputs[
"pconfig"].insert ( std::pair<std::string,double>(std::to_string(i),1.) );
496 inputs[
"pconfig"].insert ( std::pair<std::string,double>(std::to_string(i),0.) );
499 for (
int i = 0; i< 2; i ++){
501 inputs[
"econfig"].insert ( std::pair<std::string,double>(std::to_string(i),1.) );
504 inputs[
"econfig"].insert ( std::pair<std::string,double>(std::to_string(i),0.) );
508 inputs[
"ripos"].insert ( std::pair<std::string,double>(
"0", riImpactEta) );
509 inputs[
"ripos"].insert ( std::pair<std::string,double>(
"1", riImpactPhi ) );
511 return StatusCode::SUCCESS;
521 if (impactCellDDE==
nullptr){
522 ATH_MSG_WARNING(
"No cell found for this eta " << etaExtrap <<
" phi " << phiExtrap);
523 return StatusCode::FAILURE;
529 const int caloHashImpactCell=impactCellDDE->
calo_hash();
530 const double etaImpactCell=impactCellDDE->
eta();
531 const double phiImpactCell=impactCellDDE->
phi();
532 const double etaRawImpactCell=impactCellDDE->
eta_raw();
533 const double phiRawImpactCell=impactCellDDE->
phi_raw();
537 " eta=" << etaImpactCell <<
538 " phi=" << phiImpactCell <<
539 " eta raw=" << etaRawImpactCell <<
540 " phi raw=" << phiRawImpactCell );
555 sampling = theCell->caloDDE()->getSampling();
556 eta_raw = theCell->caloDDE()->eta_raw();
557 phi_raw = theCell->caloDDE()->phi_raw();
570 if ((sampling == 0) || (sampling == 1) ){
578 else if(sampling == 2) {
587 else if(sampling == 3){
599 return StatusCode::FAILURE;
604 if (etaRawImpactCell < 0){
611 return StatusCode::SUCCESS;
Scalar eta() const
pseudorapidity method
Scalar phi() const
phi method
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
This class provides an interface to generate or decode an identifier for the upper levels of the dete...
std::unique_ptr< CaloDetDescrManager > buildCaloDetDescrNoAlign(ISvcLocator *svcLocator, IMessageSvc *msgSvc)
Definition of CaloDetDescrManager.
static constexpr const char * caloMgrStaticKey
CaloPhiRange class declaration.
#define CHECK(...)
Evaluate an expression and check for errors.
bool compCellsForDNNSortMirror(const CaloCell *a, const CaloCell *b)
bool compCellsForDNNSort(const CaloCell *a, const CaloCell *b)
std::string PathResolverFindCalibFile(const std::string &logical_file_name)
constexpr int pow(int base, int exp) noexcept
Container class for CaloCell.
CaloSampling::CaloSample CaloSample
Data object for each calorimeter readout cell.
This class groups all DetDescr information related to a CaloCell.
IdentifierHash calo_hash() const
cell calo hash
float eta_raw() const
cell eta_raw
Identifier identify() const override final
cell identifier
float eta() const
cell eta
float phi() const
cell phi
float phi_raw() const
cell phi_raw
This class provides the client interface for accessing the detector description information common to...
This class initializes the Calo (LAr and Tile) offline identifiers.
const LArEM_ID * getEM_ID(void) const
static double diff(double phi1, double phi2)
simple phi1 - phi2 calculation, but result is fixed to respect range.
A manager class providing access to readout geometry information for the forward calorimeter.
const ServiceHandle< StoreGateSvc > & detStore() const
The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.
Gaudi::Property< std::string > m_screenOutputPrefix
Screen output refinement.
const ServiceHandle< StoreGateSvc > & evtStore() const
The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
BaseSimulationSvc(const std::string &name, ISvcLocator *pSvcLocator)
ServiceHandle< IChronoStatSvc > m_chrono
The timing service for general usage.
const CaloDetDescrManager * m_caloDetDescrManager
std::vector< CaloCell * > m_windowCells
virtual StatusCode initialize() override final
Athena algorithm's interface methods.
const double m_phiRawMiddleCut
virtual StatusCode simulate(ISFParticle &isp, McEventCollection *) override final
Simulation Call.
const int m_numberOfCellsForDNN
StatusCode fillWindowCells(const double etaExtrap, const double phiExtrap, const CaloDetDescrElement *&impactCellDDE)
std::string m_caloCellsOutputName
CaloCellContainer * m_theContainer
StatusCode fillNetworkInputs(const ISF::ISFParticle &isfp, NetworkInputs &inputs, double &trueEnergy)
std::string m_paramsInputArchitecture
DNNCaloSimSvc(const std::string &name, ISvcLocator *pSvcLocator)
Constructor with parameters.
std::string m_randomEngineName
double m_riImpactEtaScale
virtual StatusCode finalize() override final
framework methods
ServiceHandle< IAtRndmGenSvc > m_rndGenSvc
StatusCode initializeNetwork()
helper for initialize
std::unique_ptr< lwt::LightweightGraph > m_graph
virtual StatusCode releaseEvent() override final
Release Event chain - in case of an end-of event action is needed.
const double m_phiRawStripCut
double m_logTrueEnergyScale
CLHEP::HepRandomEngine * m_randomEngine
ToolHandleArray< ICaloCellMakerTool > m_caloCellMakerToolsSetup
const double m_etaRawBackCut
ToolHandle< IFastCaloSimCaloExtrapolation > m_FastCaloSimCaloExtrapolation
std::unique_ptr< CaloGeometryFromCaloDDM > m_caloGeo
std::map< std::string, std::map< std::string, double > > NetworkInputs
ToolHandleArray< ICaloCellMakerTool > m_caloCellMakerToolsRelease
std::map< std::string, double > NetworkOutputs
double m_logTrueEnergyMean
double m_riImpactPhiScale
virtual ~DNNCaloSimSvc() final
Destructor.
std::string m_paramsFilename
const double m_etaRawMiddleCut
virtual StatusCode setupEvent() override final
Setup Event chain - in case of a begin-of event action is needed.
The generic ISF particle definition,.
double ekin() const
Kinetic energy.
const Amg::Vector3D & momentum() const
The current momentum vector of the ISFParticle.
const Amg::Vector3D & position() const
The current position of the ISFParticle.
int pdgCode() const
PDG value.
double mass() const
mass of the particle
This defines the McEventCollection, which is really just an ObjectVector of McEvent objectsFile: Gene...
CLHEP::HepRandomEngine * randomEngine()
void set_vertex(const TLorentzVector &val)
singleton-like access to IMessageSvc via open function and helper
Eigen::Matrix< double, 3, 1 > Vector3D
IMessageSvc * getMessageSvc(bool quiet=false)
@ VIEW_ELEMENTS
this data object is a view, it does not own its elmts
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.