From be29fae293d459e5ca4f7a90c807633f248717ac Mon Sep 17 00:00:00 2001 From: John Kenneth Anders <john.kenneth.anders@cern.ch> Date: Thu, 7 Feb 2019 15:16:52 +0000 Subject: [PATCH] Merge branch '21.0-dnncalosimv1' into '21.0' 21.0 first version of DNNCaloSimSvc See merge request atlas/athena!20479 --- .../ISF_Config/python/ISF_ConfigConfigDb.py | 4 +- .../ISF/ISF_Config/python/ISF_MainConfig.py | 25 +- .../ISF_FastCaloSimServices/CMakeLists.txt | 5 +- .../python/ISF_FastCaloSimServicesConfig.py | 23 +- .../python/ISF_FastCaloSimServicesConfigDb.py | 3 +- .../src/DNNCaloSimSvc.cxx | 616 ++++++++++++++++++ .../src/DNNCaloSimSvc.h | 132 ++++ .../ISF_FastCaloSimServices_entries.cxx | 2 + .../python/ISF_SimulationSelectorsConfig.py | 6 +- .../python/ISF_SimulationSelectorsConfigDb.py | 3 +- 10 files changed, 811 insertions(+), 8 deletions(-) create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.cxx create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.h diff --git a/Simulation/ISF/ISF_Config/python/ISF_ConfigConfigDb.py b/Simulation/ISF/ISF_Config/python/ISF_ConfigConfigDb.py index 9e8f8cad9ac..8808d0d3b2f 100644 --- a/Simulation/ISF/ISF_Config/python/ISF_ConfigConfigDb.py +++ b/Simulation/ISF/ISF_Config/python/ISF_ConfigConfigDb.py @@ -1,4 +1,4 @@ -# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration """ Configuration database for ISF @@ -45,9 +45,11 @@ addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_G4HS_FastPileup", "ISF_Ker addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_ATLFASTIIF_IDOnly", "ISF_Kernel_ATLFASTIIF_IDOnly") addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_ATLFASTIIF_IDCalo", "ISF_Kernel_ATLFASTIIF_IDCalo") addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_G4FastCalo", "ISF_Kernel_G4FastCalo") +addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_G4FastCaloDNN", "ISF_Kernel_G4FastCaloDNN") addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_Fatras_newExtrapolation","ISF_Kernel_Fatras_newExtrapolation") addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_Fatras_newExtrapolation_IDOnly","ISF_Kernel_Fatras_newExtrapolation_IDOnly") addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_FastOnly", "ISF_Kernel_FastOnly") +addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_DNNOnly", "ISF_Kernel_DNNOnly") addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_MultiSimTest", "ISF_Kernel_MultiSimTest") addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_G4GammaCones" , "ISF_Kernel_G4GammaCones" ) addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_G4GammaCones_FastCalo" , "ISF_Kernel_G4GammaCones_FastCalo" ) diff --git a/Simulation/ISF/ISF_Config/python/ISF_MainConfig.py b/Simulation/ISF/ISF_Config/python/ISF_MainConfig.py index 7e70cf71d08..b3868c8ec6e 100644 --- a/Simulation/ISF/ISF_Config/python/ISF_MainConfig.py +++ b/Simulation/ISF/ISF_Config/python/ISF_MainConfig.py @@ -1,4 +1,4 @@ -# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration """ Tools configurations for ISF @@ -364,6 +364,22 @@ def getKernel_G4FastCalo(name="ISF_Kernel_G4FastCalo", **kwargs): simFlags.SimulationFlavour = "G4FastCalo" return getKernel_GenericSimulator(name, **kwargs) +############## Simulator: G4FastCaloDNN ############### +# like G4FastCalo, replacing FastCaloSimV2 by DNNCaloSim +def getKernel_G4FastCaloDNN(name="ISF_Kernel_G4FastCaloDNN", **kwargs): + kwargs.setdefault("ParticleBroker" , 'ISF_AFIIParticleBrokerSvc') + + kwargs.setdefault("BeamPipeSimulationSelectors", [ 'ISF_DefaultAFIIGeant4Selector' ] ) + kwargs.setdefault("IDSimulationSelectors" , [ 'ISF_DefaultAFIIGeant4Selector' ] ) + kwargs.setdefault("CaloSimulationSelectors" , [ 'ISF_MuonAFIIGeant4Selector', + 'ISF_EtaGreater5ParticleKillerSimSelector', + 'ISF_DefaultDNNCaloSimSelector' ] ) + kwargs.setdefault("MSSimulationSelectors" , [ 'ISF_DefaultAFIIGeant4Selector' ] ) + kwargs.setdefault("CavernSimulationSelectors" , [ 'ISF_DefaultParticleKillerSelector' ] ) + from G4AtlasApps.SimFlags import simFlags + simFlags.SimulationFlavour = "G4FastCalo" + return getKernel_GenericSimulator(name, **kwargs) + ############## Simulator: ATLFASTII ############### def getKernel_ATLFASTII(name="ISF_Kernel_ATLFASTII", **kwargs): kwargs.setdefault("ParticleBroker" , 'ISF_AFIIParticleBrokerSvc' ) @@ -529,6 +545,13 @@ def getKernel_FastOnly(name="ISF_Kernel_FastOnly", **kwargs): kwargs.setdefault("CavernSimulationSelectors" , [ 'ISF_DefaultParticleKillerSelector' ] ) return getKernel_GenericSimulatorNoG4(name, **kwargs) +############## Simulator: DNNOnly ############### +# run DNNCaloSim standalone, for faster tests +def getKernel_DNNOnly(name="ISF_Kernel_DNNOnly", **kwargs): + kwargs.setdefault("CaloSimulationSelectors" , [ 'ISF_DefaultDNNCaloSimSelector' ] ) + return getKernel_GenericSimulatorNoG4(name, **kwargs) + + ############## Simulator: G4GammaCones ############### def getKernel_G4GammaCones(name="ISF_Kernel_G4GammaCones", **kwargs): kwargs.setdefault("BeamPipeSimulationSelectors" , [ 'ISF_PhotonConeGeant4Selector' , diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/CMakeLists.txt b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/CMakeLists.txt index 8e82ab2ff1a..96b650d200b 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/CMakeLists.txt +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/CMakeLists.txt @@ -34,13 +34,14 @@ atlas_depends_on_subdirs( PUBLIC # External dependencies: find_package( CLHEP ) find_package( HepMC ) +find_package(lwtnn) # Component(s) in the package: atlas_add_component( ISF_FastCaloSimServices src/*.cxx src/components/*.cxx - INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${HEPMC_INCLUDE_DIRS} - LINK_LIBRARIES ${CLHEP_LIBRARIES} ${HEPMC_LIBRARIES} AthenaBaseComps AthenaKernel GaudiKernel IdDictParser ISF_Interfaces TrkEventPrimitives TrkExInterfaces CaloEvent StoreGateLib SGtests NavFourMom GeneratorObjects FastCaloSimLib ISF_Event ISF_FastCaloSimEvent ISF_FastCaloSimInterfaces ISF_FastCaloSimParametrizationLib PathResolver) + INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${HEPMC_INCLUDE_DIRS} ${LWTNN_INCLUDE_DIRS} + LINK_LIBRARIES ${CLHEP_LIBRARIES} ${HEPMC_LIBRARIES} ${LWTNN_LIBRARIES} AthenaBaseComps AthenaKernel GaudiKernel IdDictParser ISF_Interfaces TrkEventPrimitives TrkExInterfaces CaloEvent StoreGateLib SGtests NavFourMom GeneratorObjects FastCaloSimLib ISF_Event ISF_FastCaloSimEvent ISF_FastCaloSimInterfaces ISF_FastCaloSimParametrizationLib PathResolver) # Install files from the package: atlas_install_headers( ISF_FastCaloSimServices ) diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfig.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfig.py index 62e8ded1c78..5c476d2d9a0 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfig.py +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfig.py @@ -1,4 +1,4 @@ -# Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration """ Tools configurations for ISF_FastCaloSimServices @@ -87,3 +87,24 @@ def getFastCaloSimSvcV2(name="ISF_FastCaloSimSvcV2", **kwargs): kwargs.setdefault("RandomSvc" , simFlags.RandomSvc.get_Value() ) return CfgMgr.ISF__FastCaloSimSvcV2(name, **kwargs ) + +#### DNNCaloSim +def getDNNCaloSimSvc(name="ISF_DNNCaloSimSvc", **kwargs): + from ISF_FastCaloSimServices.ISF_FastCaloSimJobProperties import ISF_FastCaloSimFlags + + kwargs.setdefault("CaloCellsOutputName" , ISF_FastCaloSimFlags.CaloCellsName() ) + kwargs.setdefault("CaloCellMakerTools_setup" , [ 'ISF_EmptyCellBuilderTool' ] ) + kwargs.setdefault("CaloCellMakerTools_release" , [ 'ISF_CaloCellContainerFinalizerTool', + 'ISF_FastHitConvertTool' ]) #DR needed ? + kwargs.setdefault("ParamsInputFilename" , ISF_FastCaloSimFlags.ParamsInputFilename()) + kwargs.setdefault("FastCaloSimCaloExtrapolation" , 'FastCaloSimCaloExtrapolation') + + # register the FastCaloSim random number streams + from G4AtlasApps.SimFlags import simFlags + if not simFlags.RandomSeedList.checkForExistingSeed(ISF_FastCaloSimFlags.RandomStreamName()): + simFlags.RandomSeedList.addSeed( ISF_FastCaloSimFlags.RandomStreamName(), 98346412, 12461240 ) + + kwargs.setdefault("RandomStream" , ISF_FastCaloSimFlags.RandomStreamName()) + kwargs.setdefault("RandomSvc" , simFlags.RandomSvc.get_Value() ) + + return CfgMgr.ISF__DNNCaloSimSvc(name, **kwargs ) diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfigDb.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfigDb.py index 70399972113..ef876dd16d2 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfigDb.py +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfigDb.py @@ -1,4 +1,4 @@ -# Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration """ Configuration database for ISF_FastCaloSimServices @@ -35,3 +35,4 @@ addService("ISF_FastCaloSimServices.ISF_FastCaloSimServicesConfig.getFastHitConv addService("ISF_FastCaloSimServices.ISF_FastCaloSimServicesConfig.getLegacyAFIIFastCaloSimSvc", "ISF_LegacyAFIIFastCaloSimSvc") addService("ISF_FastCaloSimServices.ISF_FastCaloSimServicesConfig.getFastHitConvAlgLegacyAFIIFastCaloSimSvc", "ISF_FastHitConvAlgLegacyAFIIFastCaloSimSvc") addService("ISF_FastCaloSimServices.ISF_FastCaloSimServicesConfig.getFastCaloSimSvcV2", "ISF_FastCaloSimSvcV2") +addService("ISF_FastCaloSimServices.ISF_FastCaloSimServicesConfig.getDNNCaloSimSvc", "ISF_DNNCaloSimSvc") diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.cxx new file mode 100644 index 00000000000..30499bb8379 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.cxx @@ -0,0 +1,616 @@ +/* + Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +*/ + +/////////////////////////////////////////////////////////////////// +// DNNCaloSimSvc.cxx, (c) ATLAS Detector software // +/////////////////////////////////////////////////////////////////// + +// class header include +#include "DNNCaloSimSvc.h" + + +// FastCaloSim includes +#include "ISF_FastCaloSimEvent/TFCSParametrizationBase.h" +#include "ISF_FastCaloSimEvent/TFCSSimulationState.h" +#include "ISF_FastCaloSimEvent/TFCSTruthState.h" +#include "ISF_FastCaloSimEvent/TFCSExtrapolationState.h" +#include "ISF_FastCaloSimParametrization/CaloGeometryFromCaloDDM.h" +#include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h" + +//calculators +#include "CaloGeoHelpers/CaloPhiRange.h" + +//! +#include "AtlasDetDescr/AtlasDetectorID.h" +#include "IdDictParser/IdDictParser.h" +#include "CaloIdentifier/LArEM_ID.h" +#include "CaloIdentifier/TileID.h" +//! + +// StoreGate +#include "StoreGate/StoreGateSvc.h" +#include "StoreGate/StoreGate.h" + +#include "CaloDetDescr/CaloDetDescrElement.h" +#include "CaloDetDescr/CaloDetDescrManager.h" +#include "CaloIdentifier/CaloIdManager.h" +#include "LArReadoutGeometry/FCALDetectorManager.h" + +#include "PathResolver/PathResolver.h" + +#include "lwtnn/parse_json.hh" + +#include "TFile.h" +#include <fstream> +#include "CLHEP/Random/RandGauss.h" +#include "CLHEP/Random/RandFlat.h" +#include "CLHEP/Units/SystemOfUnits.h" + +using std::abs; +using std::atan2; + +/** Constructor **/ +ISF::DNNCaloSimSvc::DNNCaloSimSvc(const std::string& name, ISvcLocator* svc) : + BaseSimulationSvc(name, svc), + m_graph(nullptr), + m_theContainer(nullptr), + m_rndGenSvc("AtRndmGenSvc", name), + m_randomEngine(nullptr), + m_caloDetDescrManager(nullptr), + m_caloGeo(nullptr) +{ + declareProperty("ParamsInputFilename" , m_paramsFilename="DNNCaloSim/DNNCaloSim_GAN_nn_v0.json"," lwtnn json output describing the trained network"); + declareProperty("ParamsInputArchitecture" , m_paramsInputArchitecture="GANv0","tag describing additional parameters necessary for the network"); + declareProperty("CaloCellsOutputName" , m_caloCellsOutputName) ; + declareProperty("CaloCellMakerTools_setup" , m_caloCellMakerToolsSetup) ; + declareProperty("CaloCellMakerTools_release" , m_caloCellMakerToolsRelease) ; + declareProperty("RandomSvc" , m_rndGenSvc ); + declareProperty("RandomStream" , m_randomEngineName ); + declareProperty("FastCaloSimCaloExtrapolation" , m_FastCaloSimCaloExtrapolation ); +} + +ISF::DNNCaloSimSvc::~DNNCaloSimSvc() +{} + +/** framework methods */ +StatusCode ISF::DNNCaloSimSvc::initialize() +{ + ATH_MSG_INFO(m_screenOutputPrefix << "Initializing ..."); + + ATH_CHECK(m_rndGenSvc.retrieve()); + m_randomEngine = m_rndGenSvc->GetEngine( m_randomEngineName); + if(!m_randomEngine) + { + ATH_MSG_ERROR("Could not get random number engine from RandomNumberService. Abort."); + return StatusCode::FAILURE; + } + + m_caloDetDescrManager = CaloDetDescrManager::instance(); + const FCALDetectorManager * fcalManager=NULL; + ATH_CHECK(detStore()->retrieve(fcalManager)); + + const CaloIdManager* caloId_mgr = m_caloDetDescrManager->getCalo_Mgr(); + m_emID = caloId_mgr->getEM_ID(); + + m_caloGeo = std::make_unique<CaloGeometryFromCaloDDM>(); + m_caloGeo->LoadGeometryFromCaloDDM(m_caloDetDescrManager); + if(!m_caloGeo->LoadFCalChannelMapFromFCalDDM(fcalManager) )ATH_MSG_FATAL("Found inconsistency between FCal_Channel map and GEO file. Please, check if they are configured properly."); + + + // initialize DNN + if (initializeNetwork().isFailure()) + { + ATH_MSG_ERROR("Could not initialize network "); + return StatusCode::FAILURE; + + } + + + // Get FastCaloSimCaloExtrapolation + if(m_FastCaloSimCaloExtrapolation.retrieve().isFailure()) + { + ATH_MSG_ERROR("FastCaloSimCaloExtrapolation not found "); + return StatusCode::FAILURE; + } + + m_windowCells.reserve(m_numberOfCellsForDNN); + + + + return StatusCode::SUCCESS; +} + +// initialize lwtnn network +StatusCode ISF::DNNCaloSimSvc::initializeNetwork() +{ + + + + // get neural net JSON file as an std::istream object + std::string inputFile=PathResolverFindCalibFile(m_paramsFilename); + if (inputFile.empty()){ + ATH_MSG_ERROR("Could not find json file " << m_paramsFilename ); + return StatusCode::FAILURE; + } + + + + ATH_MSG_INFO("Using json file " << inputFile ); + std::ifstream input(inputFile); + // build the graph + m_graph=std::make_unique<lwt::LightweightGraph>(lwt::parse_json_graph(input) ); + if (m_graph==nullptr){ + ATH_MSG_ERROR("Could not create LightWeightGraph from " << m_paramsFilename ); + return StatusCode::FAILURE; + } + + + + // initialize all necessary constants + // FIXME eventually all these could be stored in the .json file + + ATH_MSG_INFO("Using ParamsInputArchitecture: " << m_paramsInputArchitecture ); + if (m_paramsInputArchitecture=="GANv0") // GAN then VAE etc... + { + m_GANLatentSize = 300; + m_logTrueEnergyMean = 9.70406053; + m_logTrueEnergyScale = 1.76099569; + m_riImpactEtaMean = 3.47603256e-05; + m_riImpactEtaScale = 0.00722316; + m_riImpactPhiMean = -5.42153684e-05; + m_riImpactPhiScale = 0.00708241; + } + + if (m_GANLatentSize==0){ + ATH_MSG_ERROR("m_GANLatentSize uninitialized!.ParamsInputArchitecture= " << m_paramsInputArchitecture ); + return StatusCode::FAILURE; + } + + return StatusCode::SUCCESS; +} + +/** framework methods */ +StatusCode ISF::DNNCaloSimSvc::finalize() +{ + ATH_MSG_INFO(m_screenOutputPrefix << "Finalizing ..."); + return StatusCode::SUCCESS; +} + +StatusCode ISF::DNNCaloSimSvc::setupEvent() +{ + const EventContext& ctx = Gaudi::Hive::currentContext(); + ATH_MSG_INFO(m_screenOutputPrefix << "setupEvent NEW EVENT! "); + + m_theContainer = new CaloCellContainer(SG::VIEW_ELEMENTS); + + StatusCode sc = evtStore()->record(m_theContainer, m_caloCellsOutputName); + if (sc.isFailure()) + { + ATH_MSG_FATAL( m_screenOutputPrefix << "cannot record CaloCellContainer " << m_caloCellsOutputName ); + return StatusCode::FAILURE; + } + + //FIXME why retrieve every event ? + CHECK( m_caloCellMakerToolsSetup.retrieve() ); + ATH_MSG_DEBUG( "Successfully retrieve CaloCellMakerTools: " << m_caloCellMakerToolsSetup ); + ToolHandleArray<ICaloCellMakerTool>::iterator itrTool = m_caloCellMakerToolsSetup.begin(); + ToolHandleArray<ICaloCellMakerTool>::iterator endTool = m_caloCellMakerToolsSetup.end(); + for (; itrTool != endTool; ++itrTool) + { + std::string chronoName=this->name()+"_"+ itrTool->name(); + if (m_chrono) m_chrono->chronoStart(chronoName); + StatusCode sc = (*itrTool)->process(m_theContainer, ctx); + if (m_chrono) { + m_chrono->chronoStop(chronoName); + ATH_MSG_DEBUG( m_screenOutputPrefix << "Chrono stop : delta " << m_chrono->chronoDelta (chronoName,IChronoStatSvc::USER) * CLHEP::microsecond / CLHEP::second << " second " ); + } + + if (sc.isFailure()) + { + ATH_MSG_ERROR( m_screenOutputPrefix << "Error executing tool " << itrTool->name() ); + return StatusCode::FAILURE; + } + } + + + // FIXME just for debugging + // exercise window building + // should be false by default + // careful:if enabled change the random number sequence ! + if (false){ + TFCSSimulationState testsimulstate(m_randomEngine); + const CaloDetDescrElement * testImpactCellDDE; + const int ntrial=100; + ATH_MSG_INFO ("Trial window building on " << ntrial << " dummy eta phi " ); + for (int i=0 ; i< ntrial ; i++){ + const double eta = CLHEP::RandFlat::shoot(testsimulstate.randomEngine(), 0.2, 0.25); + const double phi = CLHEP::RandFlat::shoot(testsimulstate.randomEngine(), -CLHEP::pi , CLHEP::pi); + + //randomise eta, phi + if (fillWindowCells(eta,phi,testImpactCellDDE).isFailure()){ + ATH_MSG_WARNING("Could not build trial window cells vector with eta " << eta << " phi " << phi); + } + } + ATH_MSG_INFO ("End of trial window building on " << ntrial << " dummy eta phi " ); + + } + + + + + return StatusCode::SUCCESS; +} + +StatusCode ISF::DNNCaloSimSvc::releaseEvent() +{ + const EventContext& ctx = Gaudi::Hive::currentContext(); + ATH_MSG_VERBOSE(m_screenOutputPrefix << "release Event"); + + CHECK( m_caloCellMakerToolsRelease.retrieve() ); + + //run release tools in a loop + ToolHandleArray<ICaloCellMakerTool>::iterator itrTool = m_caloCellMakerToolsRelease.begin(); + ToolHandleArray<ICaloCellMakerTool>::iterator endTool = m_caloCellMakerToolsRelease.end(); + for (; itrTool != endTool; ++itrTool) + { + ATH_MSG_VERBOSE( m_screenOutputPrefix << "Calling tool " << itrTool->name() ); + + StatusCode sc = (*itrTool)->process(m_theContainer, ctx); + + if (sc.isFailure()) + { + ATH_MSG_ERROR( m_screenOutputPrefix << "Error executing tool " << itrTool->name() ); + } + } + + return StatusCode::SUCCESS; + +} +bool compCellsForDNNSortMirror(const CaloCell* a, const CaloCell* b) +{ + static CaloPhiRange range; + CaloCell_ID::CaloSample aSampling = a->caloDDE()->getSampling(); + CaloCell_ID::CaloSample bSampling = b->caloDDE()->getSampling(); + + float aEtaRaw = a->caloDDE()->eta_raw(); + float bEtaRaw = b->caloDDE()->eta_raw(); + + float aPhiRaw = a->caloDDE()->phi_raw(); + float bPhiRaw = b->caloDDE()->phi_raw(); + + if ((aSampling) < (bSampling)) + return true; + else if ((aSampling) > (bSampling)) + return false; + // reverse sort in eta for left half of detector + if ((aEtaRaw) < (bEtaRaw)) + return false; + else if ((aEtaRaw) > (bEtaRaw)) + return true; + + if (range.diff(aPhiRaw, bPhiRaw) > 0){ + return false; + } + + return true; +} + +bool compCellsForDNNSort(const CaloCell* a, const CaloCell* b) +{ + static CaloPhiRange range; + CaloCell_ID::CaloSample aSampling = a->caloDDE()->getSampling(); + CaloCell_ID::CaloSample bSampling = b->caloDDE()->getSampling(); + + float aEtaRaw = a->caloDDE()->eta_raw(); + float bEtaRaw = b->caloDDE()->eta_raw(); + + float aPhiRaw = a->caloDDE()->phi_raw(); + float bPhiRaw = b->caloDDE()->phi_raw(); + + if ((aSampling) < (bSampling)) + return true; + else if ((aSampling) > (bSampling)) + return false; + + if ((aEtaRaw) < (bEtaRaw)) + return true; + else if ((aEtaRaw) > (bEtaRaw)) + return false; + +if (range.diff(aPhiRaw, bPhiRaw) > 0){ + return false; + } + + return true; +} + +/** Simulation Call */ +StatusCode ISF::DNNCaloSimSvc::simulate(const ISF::ISFParticle& isfp, McEventCollection*) +{ + + ATH_MSG_VERBOSE("NEW PARTICLE! DNNCaloSimSvc called with ISFParticle: " << isfp); + + + //Don't simulate particles with total energy below 10 MeV + if(isfp.ekin() < 10) { + ATH_MSG_VERBOSE("Skipping particle with Ekin: " << isfp.ekin() <<" MeV. Below the 10 MeV threshold."); + return StatusCode::SUCCESS; + } + + + //Compute all inputs to the network + NetworkInputs inputs; + double trueEnergy; + if (fillNetworkInputs(isfp,inputs,trueEnergy).isFailure()) { + ATH_MSG_WARNING("Could not initialize network "); + // bail out but do not stop the job + return StatusCode::SUCCESS; + } + + + // compute the network output values + ATH_MSG_VERBOSE("neural network input = "<<inputs); + NetworkOutputs outputs = m_graph->compute(inputs); + ATH_MSG_VERBOSE("neural network output = "<<outputs); + + + // add network output energy in the right place in the calorimeter + int itr = 0; + for ( auto & windowCell : m_windowCells) { + windowCell->addEnergy(trueEnergy * outputs[std::to_string(itr)]); + itr++; + + ATH_MSG_VERBOSE(" cell eta_raw " << windowCell->caloDDE()->eta_raw() + << " phi_raw " << windowCell->caloDDE()->phi_raw() + << " sampling " << windowCell->caloDDE()->getSampling() + << " energy " << windowCell->energy()); + } + + + + + + return StatusCode::SUCCESS; +} + + + +// compute all the necessary inputs to the network +StatusCode ISF::DNNCaloSimSvc::fillNetworkInputs(const ISF::ISFParticle& isfp, NetworkInputs & inputs, double & trueEnergy) +{ + static CaloPhiRange range; + Amg::Vector3D particle_position = isfp.position(); + Amg::Vector3D particle_direction(isfp.momentum().x(),isfp.momentum().y(),isfp.momentum().z()); + + //int barcode=isfp.barcode(); // isfp barcode, eta and phi: in case we need them + // float eta_isfp = particle_position.eta(); + // float phi_isfp = particle_position.phi(); + + TFCSTruthState truth(isfp.momentum().x(),isfp.momentum().y(),isfp.momentum().z(),sqrt(isfp.momentum().mag2()+pow(isfp.mass(),2)),isfp.pdgCode()); + truth.set_vertex(particle_position[Amg::x], particle_position[Amg::y], particle_position[Amg::z]); + + trueEnergy = isfp.ekin(); + + + TFCSExtrapolationState extrapol; + //FIXME this is extrapolating to many layers, when we only need middle layer middle surface + //FIXME could have dedicated extrapolation to save time + m_FastCaloSimCaloExtrapolation->extrapolate(extrapol,&truth); + // extrapol.Print(); + + + if (false) + { + for (int isam=0; isam< CaloCell_ID_FCS::MaxSample ; isam++){ + //enum SUBPOS { SUBPOS_MID = 0, SUBPOS_ENT = 1, SUBPOS_EXT = 2}; //MID=middle, ENT=entrance, EXT=exit of cal layer + + for (int isubpos=0; isubpos< 3 ; isubpos++){ + + ATH_MSG_VERBOSE("EXTRAPO isam=" << isam << + " isubpos=" << isubpos << + " OK=" << extrapol.OK(isam,isubpos) << + " eta=" << extrapol.eta(isam,isubpos) << + " phi=" << extrapol.phi(isam,isubpos) << + " r=" << extrapol.r(isam,isubpos) ); + + } + } + } + + //FIXME deal with endcap as well + int isam=CaloCell_ID_FCS::EMB2; + int isubpos=SUBPOS_ENT; + double etaExtrap=-999.; + double phiExtrap=-999.; + if (extrapol.eta(isam,isubpos)) { + etaExtrap=extrapol.eta(isam,isubpos); + phiExtrap=extrapol.phi(isam,isubpos); + } + + ATH_MSG_VERBOSE("Will use isam=" << isam << + " isubpos=" << isubpos << + " eta=" << etaExtrap << + " phi=" << phiExtrap ); + + + + // fill vector of cells, and DDE of middle cell + const CaloDetDescrElement * impactCellDDE; + if (fillWindowCells(etaExtrap,phiExtrap,impactCellDDE).isFailure()){ + ATH_MSG_WARNING("Could not build window cells vector "); + return StatusCode::FAILURE; + } + + + // start neural network part + // fill a map of input nodes + // this is for GAN + // most likely it should be specialised as a function of m_ParamsInputArchitecture + + double randGaussz = 0.; + + + int impactEtaIndex = m_emID->eta(impactCellDDE->identify()); + int impactPhiIndex = m_emID->phi(impactCellDDE->identify()); + double etaRawImpactCell=impactCellDDE->eta_raw(); + double phiRawImpactCell=impactCellDDE->phi_raw(); + + ATH_MSG_VERBOSE("impact eta_index " << impactEtaIndex + << " phi_index " << impactPhiIndex + << " sampling " << m_emID->sampling(impactCellDDE->identify())); + + int pconf = impactPhiIndex % 4 ; + int econf = (impactEtaIndex + 1) % 2 ; // ofset corresponds to difference in index calculated for neural net preprocessing + + double riImpactEta = (etaExtrap - etaRawImpactCell); + // mirror for left half + if (etaExtrap < 0){ + riImpactEta *= -1.; + + } + // scale & centre + riImpactEta = (riImpactEta - m_riImpactEtaMean)/m_riImpactEtaScale; + double riImpactPhi = (range.diff(phiExtrap, phiRawImpactCell) - m_riImpactPhiMean)/m_riImpactPhiScale; + + + // fill randomize latent space + TFCSSimulationState simulstate(m_randomEngine); + + //FIXME generate in one go + for (int i = 0; i< m_GANLatentSize; i ++) + { + randGaussz = CLHEP::RandGauss::shoot(simulstate.randomEngine(), 0., 1.); + inputs["Z"].insert ( std::pair<std::string,double>(std::to_string(i), randGaussz) ); + + } + + // fill preprocessed true energy + inputs["E_true"].insert ( std::pair<std::string,double>("0", (std::log(trueEnergy) - m_logTrueEnergyMean)/m_logTrueEnergyScale) ); + // fill p,e configurations multi-hot vector + for (int i = 0; i< 4; i ++) + { + if (i == pconf){ + inputs["pconfig"].insert ( std::pair<std::string,double>(std::to_string(i),1.) ); + } + else{ + inputs["pconfig"].insert ( std::pair<std::string,double>(std::to_string(i),0.) ); + } + } + for (int i = 0; i< 2; i ++){ + if (i == econf){ + inputs["econfig"].insert ( std::pair<std::string,double>(std::to_string(i),1.) ); + } + else{ + inputs["econfig"].insert ( std::pair<std::string,double>(std::to_string(i),0.) ); + } + } + // fill position of extrap particle in impact cell + inputs["ripos"].insert ( std::pair<std::string,double>("0", riImpactEta) ); + inputs["ripos"].insert ( std::pair<std::string,double>("1", riImpactPhi ) ); + + return StatusCode::SUCCESS; +} + + +StatusCode ISF::DNNCaloSimSvc::fillWindowCells(const double etaExtrap,const double phiExtrap,const CaloDetDescrElement* & impactCellDDE) { + + //now find the cell it corresponds to + //FIXME this is barrel should also look in endcap + // (note that is really looking up eta, phi, not raw eta phi + static CaloPhiRange range; + + impactCellDDE=m_caloDetDescrManager->get_element(CaloCell_ID::EMB2,etaExtrap,phiExtrap); + if (impactCellDDE==nullptr){ + ATH_MSG_WARNING("No cell found for this eta " << etaExtrap << " phi " << phiExtrap); + return StatusCode::FAILURE; + } + + + + + const int caloHashImpactCell=impactCellDDE->calo_hash(); + const double etaImpactCell=impactCellDDE->eta(); + const double phiImpactCell=impactCellDDE->phi(); + const double etaRawImpactCell=impactCellDDE->eta_raw(); + const double phiRawImpactCell=impactCellDDE->phi_raw(); + + + ATH_MSG_VERBOSE("impact cell calohash=" << caloHashImpactCell << + " eta=" << etaImpactCell << + " phi=" << phiImpactCell << + " eta raw=" << etaRawImpactCell << + " phi raw=" << phiRawImpactCell ); + + + int nSqCuts = 0; + + + // select the cells DNN will simulate + // note that m_theCellContainer has all the calorimeter cells + // this will hold the list of cells to simulate + //FIXME this can be sped up certainly + m_windowCells.clear(); + CaloCell_ID::CaloSample sampling; + float eta_raw; + float phi_raw; + for(const auto& theCell : * m_theContainer) { + sampling = theCell->caloDDE()->getSampling(); + eta_raw = theCell->caloDDE()->eta_raw(); + phi_raw = theCell->caloDDE()->phi_raw(); + if (( eta_raw < etaRawImpactCell + m_etaRawBackCut) && (eta_raw > etaRawImpactCell - m_etaRawBackCut)) { + if ((range.diff(phi_raw, phiRawImpactCell) < m_phiRawStripCut ) && (range.diff(phi_raw, phiRawImpactCell) > - m_phiRawStripCut) ) { + + } + else{ + continue; + } + } + else{ + continue; + } + + if ((sampling == 0) || (sampling == 1) ){ + if ((eta_raw < etaRawImpactCell + m_etaRawMiddleCut) && (eta_raw > etaRawImpactCell - m_etaRawMiddleCut)) { + nSqCuts ++; + // add to vector + m_windowCells.push_back(theCell); + + } + } + else if((sampling == 2)) { + if ((range.diff(phi_raw , phiRawImpactCell) < m_phiRawMiddleCut) && (range.diff(phi_raw, phiRawImpactCell) > - m_phiRawMiddleCut)) { + if ((eta_raw < etaRawImpactCell + m_etaRawMiddleCut) && (eta_raw > etaRawImpactCell - m_etaRawMiddleCut)) { + nSqCuts ++; + m_windowCells.push_back(theCell); + } + } + } + + else if(sampling == 3){ + if ((range.diff(phi_raw, phiRawImpactCell) < m_phiRawMiddleCut) && (range.diff(phi_raw , phiRawImpactCell) > - m_phiRawMiddleCut )) { + nSqCuts ++; + m_windowCells.push_back(theCell); + } + + } + } + + if (nSqCuts != m_numberOfCellsForDNN){ + ATH_MSG_WARNING("Total cells passing DNN selection is " << nSqCuts << " but should be " << m_numberOfCellsForDNN ); + // bail out, but do not stop the job + return StatusCode::FAILURE; + + } + + //sort cells within the cluster like they are fed to DNN + if (etaRawImpactCell < 0){ + std::sort(m_windowCells.begin(), m_windowCells.end(), &compCellsForDNNSortMirror); + } + else{ + std::sort(m_windowCells.begin(), m_windowCells.end(), &compCellsForDNNSort); + } + + return StatusCode::SUCCESS; + +} diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.h new file mode 100644 index 00000000000..03f75dd0ec9 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.h @@ -0,0 +1,132 @@ +/* + Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +*/ + +////////////////////////////////////////////////////////////////// +// DNNCaloSimSvc.h, (c) ATLAS Detector software +/////////////////////////////////////////////////////////////////// + +#ifndef ISF_DNNCALOSIMSVC_H +#define ISF_DNNCALOSIMSVC_H 1 + +// ISF includes +#include "ISF_Interfaces/BaseSimulationSvc.h" + +// Framework includes +#include "GaudiKernel/IChronoStatSvc.h" + +// FastCaloSim includes +#include "ISF_FastCaloSimParametrization/IFastCaloSimCaloExtrapolation.h" +#include "TrkExInterfaces/ITimedExtrapolator.h" + +#include "CaloInterface/ICaloCellMakerTool.h" + +#include "AthenaKernel/IAtRndmGenSvc.h" + +#include "AtlasDetDescr/AtlasDetectorID.h" +#include "CaloIdentifier/LArEM_ID.h" +#include "CaloIdentifier/LArHEC_ID.h" +#include "CaloIdentifier/LArFCAL_ID.h" +#include "CaloIdentifier/TileID.h" +#include "CaloDetDescr/CaloDetDescrManager.h" + +#include "lwtnn/LightweightGraph.hh" +#include "CaloEvent/CaloCellContainer.h" +#include "CLHEP/Units/PhysicalConstants.h" + +namespace CLHEP +{ + class HepRandomEngine; +} + +//forward declarations +class CaloCellContainer; +class CaloGeometryFromCaloDDM; +class TFCSParametrizationBase; + + +namespace ISF { + /** @class DNNCaloSimSvc + + @author Aishik.Ghosh -at- cern.ch, David Rousseau -at- cern.ch, + */ + + class DNNCaloSimSvc : public BaseSimulationSvc + { + public: + /** Constructor with parameters */ + DNNCaloSimSvc(const std::string& name, ISvcLocator* pSvcLocator); + + /** Destructor */ + virtual ~DNNCaloSimSvc() final; + + /** Athena algorithm's interface methods */ + virtual StatusCode initialize() override final; + virtual StatusCode finalize() override final; + + /** helper for initialize */ + StatusCode initializeNetwork(); + + /** Simulation Call */ + virtual StatusCode simulate(const ISFParticle& isp, McEventCollection*) override final; + // type of input requested by lwtnn + typedef std::map<std::string, std::map<std::string, double> > NetworkInputs ; + typedef std::map<std::string, double> NetworkOutputs; + StatusCode fillNetworkInputs(const ISF::ISFParticle& isfp, NetworkInputs & inputs, double & trueEnergy); + StatusCode fillWindowCells(const double etaExtrap,const double phiExtrap,const CaloDetDescrElement* & impactCellDDE); + + /** Setup Event chain - in case of a begin-of event action is needed */ + virtual StatusCode setupEvent() override final; + + /** Release Event chain - in case of an end-of event action is needed */ + virtual StatusCode releaseEvent() override final; + + std::string m_paramsFilename; + std::string m_paramsInputArchitecture; + + std::unique_ptr<lwt::LightweightGraph> m_graph; + + ToolHandleArray<ICaloCellMakerTool> m_caloCellMakerToolsSetup ; + ToolHandleArray<ICaloCellMakerTool> m_caloCellMakerToolsRelease ; + + ToolHandle<IFastCaloSimCaloExtrapolation> m_FastCaloSimCaloExtrapolation; + ToolHandle<Trk::ITimedExtrapolator> m_extrapolator; + + CaloCellContainer * m_theContainer; + + ServiceHandle<IAtRndmGenSvc> m_rndGenSvc; + CLHEP::HepRandomEngine* m_randomEngine; + std::string m_randomEngineName; + + + const CaloDetDescrManager* m_caloDetDescrManager; + std::unique_ptr<CaloGeometryFromCaloDDM> m_caloGeo; + const LArEM_ID* m_emID; + std::vector<CaloCell*> m_windowCells; + + // specific to architecture + // preprocessing of input + int m_GANLatentSize = 0; + double m_logTrueEnergyMean = 0.; + double m_logTrueEnergyScale = 0.; + double m_riImpactEtaMean = 0.; + double m_riImpactEtaScale = 0.; + double m_riImpactPhiMean = 0.; + double m_riImpactPhiScale = 0.; + + // building of the 266 cells cluster + const int m_numberOfCellsForDNN = 266; + const double m_middleCellWidthEta = 0.025; + const double m_middleCellWidthPhi = CLHEP::pi / std::pow(2,7); + const double m_etaRawMiddleCut = m_middleCellWidthEta * 3.5; + const double m_etaRawBackCut = m_middleCellWidthEta * 4.; + const double m_phiRawMiddleCut = m_middleCellWidthPhi * 3.5; + const double m_phiRawStripCut = m_middleCellWidthPhi * 6.0; + + + std::string m_caloCellsOutputName; + }; + +} + +#endif //> !ISF_DNNCALOSIMSVC_H diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/components/ISF_FastCaloSimServices_entries.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/components/ISF_FastCaloSimServices_entries.cxx index 703317432b7..8c8150c4f6f 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/components/ISF_FastCaloSimServices_entries.cxx +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/components/ISF_FastCaloSimServices_entries.cxx @@ -1,8 +1,10 @@ +#include "../DNNCaloSimSvc.h" #include "../FastCaloSimSvc.h" #include "../FastCaloSimSvcV2.h" #include "../FastCaloSimSvcPU.h" #include "../FastCaloTool.h" +DECLARE_COMPONENT( ISF::DNNCaloSimSvc ) DECLARE_COMPONENT( ISF::FastCaloSimSvc ) DECLARE_COMPONENT( ISF::FastCaloSimSvcV2 ) DECLARE_COMPONENT( ISF::FastCaloSimSvcPU ) diff --git a/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfig.py b/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfig.py index fe3e9f4fd19..5b9304d4198 100644 --- a/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfig.py +++ b/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfig.py @@ -1,4 +1,4 @@ -# Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration """ ISF_SimulationSelectors configurations for ISF Elmar Ritsch, 04/02/2013 @@ -63,6 +63,10 @@ def getDefaultFastCaloSimV2Selector(name="ISF_DefaultFastCaloSimV2Selector", **k kwargs.setdefault('SimulationFlavor', SimulationFlavor.FastCaloSimV2) return getDefaultSimSelector(name, **kwargs ) +def getDefaultDNNCaloSimSelector(name="ISF_DefaultDNNCaloSimSelector", **kwargs): + kwargs.setdefault("Simulator" , 'ISF_DNNCaloSimSvc') + return getDefaultSimSelector(name, **kwargs ) + def getFastHitConvAlgFastCaloSimSelector(name="ISF_FastHitConvAlgFastCaloSimSelector", **kwargs): kwargs.setdefault("Simulator" , 'ISF_FastHitConvAlgFastCaloSimSvc') kwargs.setdefault('SimulationFlavor', SimulationFlavor.FastCaloSim) diff --git a/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfigDb.py b/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfigDb.py index f64e7493ed7..d166e8912ba 100644 --- a/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfigDb.py +++ b/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfigDb.py @@ -1,4 +1,4 @@ -# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration """ Configuration database for ISF_SimulationSelectors Elmar Ritsch, 10/11/2014 @@ -36,6 +36,7 @@ addTool("ISF_SimulationSelectors.ISF_SimulationSelectorsConfig.getFastCaloSimPil addTool("ISF_SimulationSelectors.ISF_SimulationSelectorsConfig.getFastHitConvAlgFastCaloSimSelector" , "ISF_FastHitConvAlgFastCaloSimSelector" ) addTool("ISF_SimulationSelectors.ISF_SimulationSelectorsConfig.getDefaultLegacyAFIIFastCaloSimSelector" , "ISF_DefaultLegacyAFIIFastCaloSimSelector") addTool("ISF_SimulationSelectors.ISF_SimulationSelectorsConfig.getDefaultFastCaloSimV2Selector" , "ISF_DefaultFastCaloSimV2Selector" ) +addTool("ISF_SimulationSelectors.ISF_SimulationSelectorsConfig.getDefaultDNNCaloSimSelector" , "ISF_DefaultDNNCaloSimSelector" ) addTool("ISF_SimulationSelectors.ISF_SimulationSelectorsConfig.getFastHitConvAlgLegacyAFIIFastCaloSimSelector" , "ISF_FastHitConvAlgLegacyAFIIFastCaloSimSelector") addTool("ISF_SimulationSelectors.ISF_SimulationSelectorsConfig.getDefaultFatrasSelector" , "ISF_DefaultFatrasSelector" ) addTool("ISF_SimulationSelectors.ISF_SimulationSelectorsConfig.getDefaultFatrasNewExtrapolationSelector", "ISF_DefaultFatrasNewExtrapolationSelector") -- GitLab