diff --git a/Simulation/ISF/ISF_Config/python/ISF_ConfigConfigDb.py b/Simulation/ISF/ISF_Config/python/ISF_ConfigConfigDb.py index 382cb1784517c1d8d47cbe7f392bf69422417030..fbe858591038650d5b32c7a5a83db2c6e6bdcd07 100644 --- a/Simulation/ISF/ISF_Config/python/ISF_ConfigConfigDb.py +++ b/Simulation/ISF/ISF_Config/python/ISF_ConfigConfigDb.py @@ -41,6 +41,7 @@ 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_FastCaloGAN", "ISF_Kernel_FastCaloGAN") addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_G4FastCaloTest", "ISF_Kernel_G4FastCaloTest") addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_G4FastCaloDNN", "ISF_Kernel_G4FastCaloDNN") addAlgorithm("ISF_Config.ISF_MainConfig.getKernel_Fatras_newExtrapolation","ISF_Kernel_Fatras_newExtrapolation") diff --git a/Simulation/ISF/ISF_Config/python/ISF_MainConfig.py b/Simulation/ISF/ISF_Config/python/ISF_MainConfig.py index 94b1e8f7d0d133e70f3d5af20ed3c60b386ccfdd..ae341c203a97aa9e52e4730b861b002fb6778829 100644 --- a/Simulation/ISF/ISF_Config/python/ISF_MainConfig.py +++ b/Simulation/ISF/ISF_Config/python/ISF_MainConfig.py @@ -352,6 +352,31 @@ def getKernel_G4FastCaloDNN(name="ISF_Kernel_G4FastCaloDNN", **kwargs): simFlags.SimulationFlavour = "G4FastCaloDNN" return getKernel_GenericSimulator(name, **kwargs) +############## Simulator: FastCaloGAN ############### +def getKernel_FastCaloGAN(name="ISF_Kernel_FastCaloGAN", **kwargs): + kwargs.setdefault("ParticleBroker" , 'ISF_AFIIParticleBrokerSvc') + + kwargs.setdefault("BeamPipeSimulationSelectors", [ 'ISF_DefaultAFIIGeant4Selector' ] ) + kwargs.setdefault("IDSimulationSelectors" , [ 'ISF_DefaultAFIIGeant4Selector' ] ) + kwargs.setdefault("CaloSimulationSelectors" , [ 'ISF_MuonAFIIGeant4Selector', + 'ISF_EtaGreater5ParticleKillerSimSelector', + 'ISF_ElectronAFIIGeant4Selector', + 'ISF_PhotonAFIIGeant4Selector', + 'ISF_PionAFIIGeant4Selector', + 'ISF_ProtonAFIIGeant4Selector', + 'ISF_ChargedKaonAFIIGeant4Selector', + 'ISF_KLongAFIIGeant4Selector', + 'ISF_DefaultDNNCaloSimSelector' ] ) + kwargs.setdefault("MSSimulationSelectors" , [ 'ISF_DefaultAFIIGeant4Selector' ] ) + kwargs.setdefault("CavernSimulationSelectors" , [ 'ISF_DefaultParticleKillerSelector' ] ) + from G4AtlasApps.SimFlags import simFlags + simFlags.SimulationFlavour = "G4FastCaloGAN" + return getKernel_G4FastCalo(name, **kwargs) + + + + + ############## Simulator: ATLFASTII ############### def getKernel_ATLFASTII(name="ISF_Kernel_ATLFASTII", **kwargs): kwargs.setdefault("ParticleBroker" , 'ISF_AFIIParticleBrokerSvc' ) diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/CMakeLists.txt b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/CMakeLists.txt index 96b650d200bece70dc7ce10ad249344c5ea76d18..965a9b80dd28cfb197b18d775dd86c2df36c52fe 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/CMakeLists.txt +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/CMakeLists.txt @@ -2,6 +2,8 @@ # Package: ISF_FastCaloSimServices ################################################################################ +set(ATHENA_FCS_PATH "${ATHENA_PATH}/Simulation/ISF/ISF_FastCaloSim") + # Declare the package name: atlas_subdir( ISF_FastCaloSimServices ) @@ -35,15 +37,21 @@ atlas_depends_on_subdirs( PUBLIC find_package( CLHEP ) find_package( HepMC ) find_package(lwtnn) +find_package( LibXml2 REQUIRED) + +include_directories( + ${LIBXML2_INCLUDE_DIR} + ) # Component(s) in the package: atlas_add_component( ISF_FastCaloSimServices src/*.cxx src/components/*.cxx - 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) + INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${HEPMC_INCLUDE_DIRS} ${LWTNN_INCLUDE_DIRS} ${LIBXML2_INCLUDE_DIR} + LINK_LIBRARIES ${CLHEP_LIBRARIES} ${HEPMC_LIBRARIES} ${LWTNN_LIBRARIES} ${LIBXML2_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 ) +atlas_install_headers( ISF_FastCaloSimParametrization ) atlas_install_python_modules( python/*.py ) diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimJobProperties.py b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimJobProperties.py index 280c5c4fd0a53908283624ee41544f837bc88d0f..721236d2257f866784b9320ef316499927d17114 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimJobProperties.py +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimJobProperties.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 ## @file ISF_FastCaloSimJobProperties.py ## @purpose Python module to hold common flags to configure JobOptions @@ -56,11 +56,29 @@ class DoRandomFluctuations(JobProperty): allowedTypes = ['bool'] StoredValue = False -class ParamsInputFilename(JobProperty): - """ Filename of the input parametrizations file. """ +class DNNCaloSimInputFilename(JobProperty): + """ Filename of the input DNNCaloSim json file. """ statusOn = True allowedTypes = ['str'] - StoredValue = 'FastCaloSim/MC16/TFCSparam_v011.root' + StoredValue = 'DNNCaloSim/DNNCaloSim_GAN_nn_v1.json' + +class DNNInputArchitecture(JobProperty): + """ Architecture of the GAN. """ + statusOn = True + allowedTypes = ['str'] + StoredValue = 'FastCaloGAN' + +class FastCaloGANInputFolderName(JobProperty): + """ Folder with all the FastCaloGAN trained GAN. """ + statusOn = True + allowedTypes = ['str'] + StoredValue = '.' + +class RunSingleGAN(JobProperty): + """ Allow to run Service if not all GANs are available. """ + statusOn = True + allowedTypes = ['bool'] + StoredValue = True ##----------------------------------------------------------------------------- ## 2nd step @@ -84,7 +102,10 @@ jobproperties.ISF_FastCaloSimJobProperties.add_JobProperty( RandomStreamName jobproperties.ISF_FastCaloSimJobProperties.add_JobProperty( CaloCellsName ) jobproperties.ISF_FastCaloSimJobProperties.add_JobProperty( FastShowerInputCollection ) jobproperties.ISF_FastCaloSimJobProperties.add_JobProperty( DoRandomFluctuations ) -jobproperties.ISF_FastCaloSimJobProperties.add_JobProperty( ParamsInputFilename ) +jobproperties.ISF_FastCaloSimJobProperties.add_JobProperty( DNNCaloSimInputFilename ) +jobproperties.ISF_FastCaloSimJobProperties.add_JobProperty( DNNInputArchitecture ) +jobproperties.ISF_FastCaloSimJobProperties.add_JobProperty( FastCaloGANInputFolderName ) +jobproperties.ISF_FastCaloSimJobProperties.add_JobProperty( RunSingleGAN ) ##----------------------------------------------------------------------------- ## 5th step 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 4c1e7d204cabc3b39528278cfe310c77734aa7f7..5bca023f9679f2f1cdef73859cf599da70688b4d 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfig.py +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/python/ISF_FastCaloSimServicesConfig.py @@ -90,7 +90,6 @@ def getFastHitConvAlgLegacyAFIIFastCaloSimSvc(name="ISF_FastHitConvAlgLegacyAFII #### FastCaloSimV2 def getFastCaloSimV2ParamSvc(name="ISF_FastCaloSimV2ParamSvc", **kwargs): from ISF_FastCaloSimServices.ISF_FastCaloSimJobProperties import ISF_FastCaloSimFlags - kwargs.setdefault("ParamsInputFilename" , ISF_FastCaloSimFlags.ParamsInputFilename()) kwargs.setdefault("ParamsInputObject" , 'SelPDGID') return CfgMgr.ISF__FastCaloSimV2ParamSvc(name, **kwargs ) @@ -125,7 +124,10 @@ def getDNNCaloSimSvc(name="ISF_DNNCaloSimSvc", **kwargs): kwargs.setdefault("CaloCellMakerTools_setup" , [ 'ISF_EmptyCellBuilderTool' ] ) kwargs.setdefault("CaloCellMakerTools_release" , [ 'ISF_CaloCellContainerFinalizerTool', 'ISF_FastHitConvertTool' ]) #DR needed ? - kwargs.setdefault("ParamsInputFilename" , ISF_FastCaloSimFlags.ParamsInputFilename()) + kwargs.setdefault("DNNCaloSimInputFilename" , ISF_FastCaloSimFlags.DNNCaloSimInputFilename()) + kwargs.setdefault("DNNInputArchitecture" , ISF_FastCaloSimFlags.DNNInputArchitecture()) + kwargs.setdefault("FastCaloGANInputFolderName" , ISF_FastCaloSimFlags.FastCaloGANInputFolderName()) + kwargs.setdefault("RunSingleGAN" , ISF_FastCaloSimFlags.RunSingleGAN()) kwargs.setdefault("FastCaloSimCaloExtrapolation" , 'FastCaloSimCaloExtrapolation') # register the FastCaloSim random number streams diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.cxx index 75bbc071cc52e9052d631668c11dd30a90e242f1..620c33cf856c595edbb47e785acae103cfa91b78 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.cxx +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.cxx @@ -46,6 +46,18 @@ #include "CLHEP/Random/RandGauss.h" #include "CLHEP/Random/RandFlat.h" +// XML reader +#include +#include +#include +#include +#include +#include + +#include "HepPDT/ParticleData.hh" +#include "HepPDT/ParticleDataTable.hh" + + using std::abs; using std::atan2; @@ -59,14 +71,16 @@ ISF::DNNCaloSimSvc::DNNCaloSimSvc(const std::string& name, ISvcLocator* svc) : 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 ); + declareProperty("DNNCaloSimInputFilename" , m_DNNCaloSimInputFilename="FastCaloGan_65GeV_photon.json","The lwtnn json describing the trained network"); // default value overwritten by python + declareProperty("DNNInputArchitecture" , m_DNNInputArchitecture="NoGan","tag describing additional parameters necessary for the network"); + declareProperty("FastCaloGANInputFolderName" , m_FastCaloGANInputFolderName=".","The location of all trained networks"); + declareProperty("RunSingleGAN" , m_RunSingleGAN=false,"If true will try to run even if some GAN are not loaded."); + 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() @@ -85,17 +99,14 @@ StatusCode ISF::DNNCaloSimSvc::initialize() return StatusCode::FAILURE; } + //Detector geometry, cell based GAN 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(); - 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."); - + m_geo = new CaloGeometryFromCaloDDM(); + m_geo->LoadGeometryFromCaloDDM(m_caloDetDescrManager); + if(!m_geo->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()) @@ -124,45 +135,81 @@ StatusCode ISF::DNNCaloSimSvc::initialize() 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::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; + ATH_MSG_INFO("Using DNNInputArchitecture: " << m_DNNInputArchitecture ); + ATH_MSG_INFO(" FastCaloGANInputFolderName: " << m_FastCaloGANInputFolderName ); + ATH_MSG_INFO(" RunSingleGAN: " << m_RunSingleGAN ); + if (m_DNNInputArchitecture=="GANv1") { + // get neural net JSON file as an std::istream object + std::string inputFile=PathResolverFindCalibFile(m_DNNCaloSimInputFilename); + if (inputFile.empty()){ + ATH_MSG_ERROR("Could not find json file " << m_DNNCaloSimInputFilename ); + return StatusCode::FAILURE; + } + + ATH_MSG_INFO("Using json file " << inputFile ); + std::ifstream input(inputFile); + // build the graph + m_graph=std::make_unique(lwt::parse_json_graph(input) ); + if (m_graph==nullptr){ + ATH_MSG_ERROR("Could not create LightWeightGraph from " << m_DNNCaloSimInputFilename ); + return StatusCode::FAILURE; } + 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; + } + else if (m_DNNInputArchitecture=="FastCaloGAN") { + // get neural net JSON file as an std::istream object + int pids[3] = {11,22,211}; + for (int i= 0; i < 3; i++ ){ + int pid = pids[i]; + for (int j = 0; j < 100; j++){ + int etaMin = j * 5; + int etaMax = etaMin + 5; + + std::string inputFileName = m_FastCaloGANInputFolderName + "/neural_net_" + std::to_string(pid) +"_eta_" + std::to_string(etaMin) +"_" + std::to_string(etaMax) +".json"; + //std::string inputFileName = m_FastCaloGANInputFolderName + "/neural_net_pions_eta_25_30.json"; + std::string inputFile=PathResolverFindCalibFile(inputFileName); + if (inputFile.empty()){ + if (!m_RunSingleGAN) { + ATH_MSG_ERROR("Could not find json file " << inputFile ); + return StatusCode::FAILURE; + } + else { + ATH_MSG_WARNING("Could not find json file " << inputFile << ", will try to load other GANs in case only a partial detector is run" ); + } + } + else { + ATH_MSG_INFO("For pid: " << pid <<" and eta " << etaMin <<"-" << etaMax <<", loading json file " << inputFile ); + std::ifstream input(inputFile); + // build the graph + auto graph=std::make_unique(lwt::parse_json_graph(input) ); + if (graph==nullptr){ + ATH_MSG_ERROR("Could not create LightWeightGraph from " << inputFile ); + return StatusCode::FAILURE; + } + + m_graphs[i][j] = std::move(graph); //Fill by position so that it can be easily retrieved + } + } + } + + m_GANLatentSize = 50; + + //Get all Binning histograms to store in memory + GetAllBinning(); + } + if (m_GANLatentSize==0){ - ATH_MSG_ERROR("m_GANLatentSize uninitialized!.ParamsInputArchitecture= " << m_paramsInputArchitecture ); + ATH_MSG_ERROR("m_GANLatentSize uninitialized!.DNNInputArchitecture= " << m_DNNInputArchitecture ); return StatusCode::FAILURE; } @@ -226,7 +273,7 @@ StatusCode ISF::DNNCaloSimSvc::setupEvent() const double phi = CLHEP::RandFlat::shoot(testsimulstate.randomEngine(), -CLHEP::pi , CLHEP::pi); //randomise eta, phi - if (fillWindowCells(eta,phi,testImpactCellDDE).isFailure()){ + if (fillWindowCells(eta,phi,testImpactCellDDE,(CaloCell_ID::CaloSample) CaloCell_ID_FCS::EMB2).isFailure()){ ATH_MSG_WARNING("Could not build trial window cells vector with eta " << eta << " phi " << phi); } } @@ -326,52 +373,404 @@ if (range.diff(aPhiRaw, bPhiRaw) > 0){ /** Simulation Call */ StatusCode ISF::DNNCaloSimSvc::simulate(const ISF::ISFParticle& isfp) { - 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; + if (m_DNNInputArchitecture=="GANv0") // GAN then VAE etc... + { + 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 = "<compute(inputs); + ATH_MSG_VERBOSE("neural network output = "<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()); + } } + else if (m_DNNInputArchitecture=="FastCaloGAN") + { + ATH_MSG_VERBOSE("Get Inputs"); + if (fillFastCaloGanNetworkInputs(isfp,inputs,trueEnergy).isFailure()) { + ATH_MSG_WARNING("Could not initialize network "); + // bail out but do not stop the job + return StatusCode::SUCCESS; + } + ATH_MSG_VERBOSE("Fill Energies"); + if (fillEnergy(isfp,inputs).isFailure()){ + ATH_MSG_WARNING("Could not fill energies "); + // bail out but do not stop the job + return StatusCode::SUCCESS; + } + } - // compute the network output values - ATH_MSG_VERBOSE("neural network input = "<compute(inputs); - ATH_MSG_VERBOSE("neural network output = "<addEnergy(trueEnergy * outputs[std::to_string(itr)]); - itr++; + trueEnergy = isfp.momentum().mag(); + double randUniformZ = 0.; + //FIXME: make dependency on input particle eta, pdgid and energy + for (int i = 0; i< m_GANLatentSize; i ++) + { + randUniformZ = CLHEP::RandFlat::shoot(simulstate.randomEngine(), -1., 1.); + inputs["node_0"].insert ( std::pair(std::to_string(i), randUniformZ) ); + } + + //std::cout << "Check label: " <("0", trueEnergy/(std::pow(2,22))) ); + return StatusCode::SUCCESS; +} + +StatusCode ISF::DNNCaloSimSvc::fillEnergy(const ISF::ISFParticle& isfp, NetworkInputs inputs) +{ + Amg::Vector3D particle_position = isfp.position(); + + 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]); - ATH_MSG_VERBOSE(" cell eta_raw " << windowCell->caloDDE()->eta_raw() - << " phi_raw " << windowCell->caloDDE()->phi_raw() - << " sampling " << windowCell->caloDDE()->getSampling() - << " energy " << windowCell->energy()); + 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); + + //The position will be calculated at the middle of each layer but the eta of the incoming particle is calculated at the entrace of the first layer + int isubpos=SUBPOS_MID; + double etaExtrap=extrapol.CaloSurface_eta(); + double phiExtrap=-999.; + double rExtrapol=-999.; + double zExtrapol=-999.; + + const int pdgId = truth.pdgid(); + const double charge = HepPDT::ParticleID(pdgId).charge(); + + if (extrapol.IDCaloBoundary_eta() < -900) { + ATH_MSG_WARNING("Extrapolator failed, skipping particle which has caloSurface eta of "<< etaExtrap); + + return StatusCode::SUCCESS; } + if(etaExtrap < -900){ + ATH_MSG_WARNING("Extrapolator failed, skipping particle"); + return StatusCode::SUCCESS; + } + if(etaExtrap < -5 || etaExtrap > 5 ){ + ATH_MSG_WARNING("Extrapolator ended up outside the detector, should not be happening as Selector for 5eta should catch this, only solution is to skipping particle"); + return StatusCode::SUCCESS; + } + ATH_MSG_VERBOSE("Momentum " << isfp.momentum().mag() <<" pdgId " << pdgId << " eta " << etaExtrap); + + int pid = GetPidIndex(pdgId); + + int eta = (fabs(etaExtrap) * 100) / 5; + + // compute the network output values + //ATH_MSG_VERBOSE("neural network input = "<compute(inputs); + //ATH_MSG_VERBOSE("neural network output = "<GetNbinsX(); + //If only one bin in r means layer is empty, no value should be added + if (xBinNum == 1) { + ATH_MSG_VERBOSE(" Layer "<< layer + << " has only one bin in r, this means is it not used, skipping (this is needed to keep correct syncronisation of voxel and layers)"); + //delete h; + continue; + } + int yBinNum = h->GetNbinsY(); + for (int iy = 1; iy <= yBinNum; ++iy){ + for (int ix = 1; ix <= xBinNum; ++ix){ + double energyInVoxel = outputs["out_" + std::to_string(vox)]; + ATH_MSG_VERBOSE(" Vox "<< vox + << " energy " << energyInVoxel + << " binx " << ix + << " biny " << iy); + + if (energyInVoxel == 0){ + vox++; + continue; + } + + TAxis* x = (TAxis*)h->GetXaxis(); + nHitsR = x->GetBinUpEdge(ix) - x->GetBinLowEdge(ix); + if (yBinNum == 1){ + //nbins in alpha depend on circumference lenght + double r = x->GetBinUpEdge(ix); + nHitsAlpha = ceil(2 * TMath::Pi() * r / binResolution); + } + else{ + //d = 2•r•sin (a/2r) this distance at the upper r must be 1mm for layer 1 or 5, 5mm otherwise. + TAxis* y = (TAxis*)h->GetYaxis(); + double angle = y->GetBinUpEdge(iy) - y->GetBinLowEdge(iy); + double r = x->GetBinUpEdge(ix); + double d = 2 * r * sin(angle/2*r); + nHitsAlpha = ceil(d/binResolution); + } + + nHitsAlpha = std::min(10,std::max(1,nHitsAlpha)); + nHitsR = std::min(10,std::max(1,nHitsR)); + + for (int ir = 0; ir < nHitsR; ++ir){ + TAxis* x = (TAxis*)h->GetXaxis(); + double r = x->GetBinLowEdge(ix) + x->GetBinWidth(ix)/(nHitsR+1)*ir; + + + for (int ialpha = 0; ialpha < nHitsAlpha; ++ialpha){ + double alpha; + if (yBinNum == 1){ + TFCSSimulationState testsimulstate(m_randomEngine); + alpha = CLHEP::RandFlat::shoot(testsimulstate.randomEngine(), -CLHEP::pi , CLHEP::pi); + } + else { + TAxis* y = (TAxis*)h->GetYaxis(); + alpha = y->GetBinLowEdge(iy) + y->GetBinWidth(iy)/(nHitsAlpha+1)*ialpha; + } + + const CaloDetDescrElement* CellDDE; + + if (layer <=20){ + float delta_eta_mm = r * cos(alpha); + float delta_phi_mm = r * sin(alpha); + + // Particles with negative eta are expected to have the same shape as those with positive eta after transformation: delta_eta --> -delta_eta + if(center_eta<0.)delta_eta_mm = -delta_eta_mm; + // Particle with negative charge are expected to have the same shape as positively charged particles after transformation: delta_phi --> -delta_phi + if(charge < 0.) delta_phi_mm = -delta_phi_mm; + + const float delta_eta = delta_eta_mm / eta_jakobi / dist000; + const float delta_phi = delta_phi_mm / center_r; + + const float cell_eta = center_eta + delta_eta; + const float cell_phi = center_phi + delta_phi; + + ATH_MSG_VERBOSE(" Hit eta " << cell_eta + << " phi " << cell_phi + << " layer " << layer + << " Sampling " << static_cast(layer)); + + //if (layer == 12) continue; + CellDDE=m_geo->getDDE(layer,cell_eta,cell_phi); + //const CaloDetDescrElement* CellDDE=m_caloDetDescrManager->get_element(static_cast(layer),cell_eta,cell_phi); + } + else { //FCAL is in (x,y,z) + const float hit_r = r*cos(alpha) + center_r; + float delta_phi = r*sin(alpha)/center_r; + // Particle with negative charge are expected to have the same shape as positively charged particles after transformation: delta_phi --> -delta_phi + if(charge < 0.) delta_phi = -delta_phi; + const float hit_phi= delta_phi + center_phi; + const float x = hit_r*cos(hit_phi); + const float y = hit_r*sin(hit_phi); + const float z = center_z; + + CellDDE=m_geo->getFCalDDE(layer,x,y,z); + } + + if (CellDDE==nullptr){ + ATH_MSG_WARNING("No cell found for layer: " << layer << " r: " << r <<" and alpha:"<calo_hash(); + const int cellIdentifyHash=CellDDE->identifyHash(); + const double etaCell=CellDDE->eta(); + const double phiCell=CellDDE->phi(); + const double etaRawCell=CellDDE->eta_raw(); + const double phiRawCell=CellDDE->phi_raw(); + double depositedEnergy = EKin*energyInVoxel/(nHitsAlpha*nHitsR); + int hitNumber = ir*nHitsAlpha+ialpha; + + ATH_MSG_VERBOSE(" Vox " << vox << " calohash=" << caloHashCell << + " cell identifyHash=" << cellIdentifyHash << + " eta=" << etaCell << + " phi=" << phiCell << + " eta raw=" << etaRawCell << + " phi raw=" << phiRawCell << + " deposited Energy=" << depositedEnergy << + " total vox Energy=" << EKin*outputs["out_" + std::to_string(vox)] << + " hit Number=" << hitNumber ); + + CaloCell* theCell = (CaloCell*) m_theContainer->findCell(caloHashCell); + theCell->addEnergy(depositedEnergy); + ATH_MSG_VERBOSE(" Vox " << vox << + " cell identifyHash=" << cellIdentifyHash << + " total deposited Energy so far=" << theCell->energy() << + " Cell eta=" << theCell->eta()); + } + } + vox++; + } + } + + ATH_MSG_VERBOSE(" Number of voxels " << vox << + " for pid " << pdgId << + " and eta " << etaExtrap); + + //delete h; + ATH_MSG_VERBOSE("Done layer "<children; nodeRoot != NULL; nodeRoot = nodeRoot->next) { + if (xmlStrEqual( nodeRoot->name, BAD_CAST "Bins" )) { + for( xmlNodePtr nodeBin = nodeRoot->children; nodeBin != NULL; nodeBin = nodeBin->next ) { + if (xmlStrEqual( nodeBin->name, BAD_CAST "Bin" )) { + int nodePid = atof( (const char*) xmlGetProp( nodeBin, BAD_CAST "pid" ) ); + //int nodeEtaMin = atof( (const char*) xmlGetProp( nodeBin, BAD_CAST "etaMin" ) ); + int nodeEtaMax = atof( (const char*) xmlGetProp( nodeBin, BAD_CAST "etaMax" ) ); + + Binning binsInLayer; + int pidIndex = GetPidIndex(nodePid); + + for( xmlNodePtr nodeLayer = nodeBin->children; nodeLayer != NULL; nodeLayer = nodeLayer->next ) { + if( xmlStrEqual( nodeLayer->name, BAD_CAST "Layer" ) ) { + std::vector edges; + std::string s( (const char*)xmlGetProp( nodeLayer, BAD_CAST "r_edges" ) ); + + std::istringstream ss(s); + std::string token; + + while(std::getline(ss, token, ',')) { + edges.push_back(atof( token.c_str() )); + } + + int binsInAlpha = atof( (const char*) xmlGetProp( nodeLayer, BAD_CAST "n_bin_alpha" ) ); + int layer = atof( (const char*) xmlGetProp( nodeLayer, BAD_CAST "id" ) ); + + ATH_MSG_INFO("Layer: "<get_element(CaloCell_ID::EMB2,etaExtrap,phiExtrap); + impactCellDDE=m_caloDetDescrManager->get_element(isam,etaExtrap,phiExtrap); if (impactCellDDE==nullptr){ ATH_MSG_WARNING("No cell found for this eta " << etaExtrap << " phi " << phiExtrap); return StatusCode::FAILURE; diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.h index c28acd50261da18cbd8d088ba084f6bd4316839f..10ad1f5fa96e11cc704c77f69f820cd442675c16 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.h +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/DNNCaloSimSvc.h @@ -33,6 +33,8 @@ #include "lwtnn/LightweightGraph.hh" #include "CaloEvent/CaloCellContainer.h" +#include "TH2D.h" + namespace CLHEP { @@ -72,8 +74,15 @@ namespace ISF { // type of input requested by lwtnn typedef std::map > NetworkInputs ; typedef std::map NetworkOutputs; + typedef std::map Binning; StatusCode fillNetworkInputs(const ISF::ISFParticle& isfp, NetworkInputs & inputs, double & trueEnergy); - StatusCode fillWindowCells(const double etaExtrap,const double phiExtrap,const CaloDetDescrElement* & impactCellDDE); + StatusCode fillWindowCells(const double etaExtrap,const double phiExtrap,const CaloDetDescrElement* & impactCellDDE,CaloCell_ID::CaloSample isam); + + /** For FastCaloGAN */ + StatusCode fillFastCaloGanNetworkInputs(const ISF::ISFParticle& isfp, NetworkInputs & inputs, double & trueEnergy); + StatusCode fillEnergy(const ISF::ISFParticle& isfp, NetworkInputs inputs); + Binning GetBinningForParticleAndEta(int pid, float eta); + void GetAllBinning(); /** Setup Event chain - in case of a begin-of event action is needed */ virtual StatusCode setupEvent() override final; @@ -81,10 +90,18 @@ namespace ISF { /** 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; + int GetPidIndex(int pid); + + std::string m_DNNCaloSimInputFilename; + std::string m_DNNInputArchitecture; + std::string m_FastCaloGANInputFolderName; + bool m_RunSingleGAN; std::unique_ptr m_graph; + std::unique_ptr m_graphs[3][100]; + + std::vector AllBinning[3]; + std::vector EtaMaxList[3]; ToolHandleArray m_caloCellMakerToolsSetup ; ToolHandleArray m_caloCellMakerToolsRelease ; @@ -104,6 +121,8 @@ namespace ISF { const LArEM_ID* m_emID; std::vector m_windowCells; + CaloGeometryFromCaloDDM* m_geo; //! do not persistify + // specific to architecture // preprocessing of input int m_GANLatentSize = 0; diff --git a/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfig.py b/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfig.py index 308544ed0cc17b90e973b395b60d47c73eaa457e..d4e7223abd74f7f489cb89de8e55eadaa241f15d 100644 --- a/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfig.py +++ b/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfig.py @@ -217,6 +217,24 @@ def getEtaGreater5PileupParticleKillerSimSelector(name="ISF_EtaGreater5PileupPar kwargs.setdefault('InvertCuts' , True ) return CfgMgr.ISF__KinematicPileupSimSelector(name, **kwargs) +def getProtonAFIIGeant4Selector(name="ISF_ProtonAFIIGeant4Selector", **kwargs): + kwargs.setdefault('MaxMom' , 750) + kwargs.setdefault('ParticlePDG' , 2112) + kwargs.setdefault('Simulator' , 'ISF_AFIIGeant4SimSvc') + return CfgMgr.ISF__KinematicSimSelector(name, **kwargs) + +def getElectronAFIIGeant4Selector(name="ISF_ElectronAFIIGeant4Selector", **kwargs): + kwargs.setdefault('MaxMom' , 200) + kwargs.setdefault('ParticlePDG' , 11) + kwargs.setdefault('Simulator' , 'ISF_AFIIGeant4SimSvc') + return CfgMgr.ISF__KinematicSimSelector(name, **kwargs) + +def getPhotonAFIIGeant4Selector(name="ISF_PhotonnAFIIGeant4Selector", **kwargs): + kwargs.setdefault('MaxMom' , 200) + kwargs.setdefault('ParticlePDG' , 22) + kwargs.setdefault('Simulator' , 'ISF_AFIIGeant4SimSvc') + return CfgMgr.ISF__KinematicSimSelector(name, **kwargs) + ### ConeSimSelector Configurations def getPhotonConeSelector(name="ISF_PhotonConeSelector", **kwargs): diff --git a/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfigDb.py b/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfigDb.py index 441640d7e994ab7d5c837ac75cda608389542337..12e2877634b4b94940fc6d41d7e3a6b8a8bb31c0 100644 --- a/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfigDb.py +++ b/Simulation/ISF/ISF_SimulationSelectors/python/ISF_SimulationSelectorsConfigDb.py @@ -61,3 +61,7 @@ addTool("ISF_SimulationSelectors.ISF_SimulationSelectorsConfig.getFatrasPileupSe addTool("ISF_SimulationSelectors.ISF_SimulationSelectorsConfig.getFatrasPileupSelector_noHits" , "ISF_FatrasPileupSelector_noHits" ) addTool("ISF_SimulationSelectors.ISF_SimulationSelectorsConfig.getFatrasRandomSelector" , "ISF_FatrasRandomSelector" ) addTool("ISF_SimulationSelectors.ISF_SimulationSelectorsConfig.getDefaultParametricSimulationSelector" , "ISF_DefaultParametricSimulationSelector" ) + +addTool("ISF_SimulationSelectors.ISF_SimulationSelectorsConfig.getElectronAFIIGeant4Selector" , "ISF_ElectronAFIIGeant4Selector" ) +addTool("ISF_SimulationSelectors.ISF_SimulationSelectorsConfig.getPhotonAFIIGeant4Selector" , "ISF_PhotonAFIIGeant4Selector" ) +