From 55886231de7ae606c8403cedff6f10ff11b1b694 Mon Sep 17 00:00:00 2001 From: John Derek Chapman <chapman@hep.phy.cam.ac.uk> Date: Fri, 22 Jul 2016 10:41:28 +0200 Subject: [PATCH] CMakeLists.txt - fixes so that all required symbols are available to client packages. ATLASSIM-3004. Tagging as ISF_FastCaloSimEvent-00-02-01 (ISF_FastCaloSimEvent-00-02-01) * CMakeLists.txt - fixes so that all required symbols are available to client packages. ATLASSIM-3004 * cmt/requirements - whitespace clean-up. * Tagging as ISF_FastCaloSimEvent-00-02-01 2015-06-01 Elmar Ritsch <Elmar.Ritsch@cern.ch> * Move following classes from ISF_FastCaloSimParametrization package into this package and make them compile-able with cmt (outside of ROOT). ATLASSIM-2882: FastCaloSim_CaloCell_ID IntArray TFCS1DFunction TFCSEnergyParametrization TFCSExtrapolationState TFCSParametrization TFCSParametrizationBase TFCSPCAEnergyParametrization TFCSSimulationState TFCSTruthState ... (Long ChangeLog diff - truncated) Former-commit-id: 0217f1af3f42617bd91e18e067552df3ed982d97 --- .../ISF_FastCaloSimEvent/CMakeLists.txt | 18 +- .../FastCaloSim_CaloCell_ID.h | 30 ++ .../ISF_FastCaloSimEvent/IntArray.h | 26 ++ .../ISF_FastCaloSimEvent/TFCS1DFunction.h | 42 ++ .../TFCSEnergyParametrization.h | 29 ++ .../TFCSExtrapolationState.h | 65 +++ .../TFCSPCAEnergyParametrization.h | 56 +++ .../TFCSParametrization.h | 49 +++ .../TFCSParametrizationBase.h | 49 +++ .../TFCSSimulationState.h | 43 ++ .../ISF_FastCaloSimEvent/TFCSTruthState.h | 33 ++ .../ISF_FastCaloSimEvent/cmt/requirements | 32 +- .../ISF_FastCaloSimEvent/src/IntArray.cxx | 31 ++ .../src/TFCS1DFunction.cxx | 371 ++++++++++++++++++ .../src/TFCSEnergyParametrization.cxx | 21 + .../src/TFCSExtrapolationState.cxx | 42 ++ .../src/TFCSPCAEnergyParametrization.cxx | 352 +++++++++++++++++ .../src/TFCSParametrization.cxx | 67 ++++ .../src/TFCSParametrizationBase.cxx | 39 ++ .../src/TFCSSimulationState.cxx | 35 ++ .../src/TFCSTruthState.cxx | 25 ++ 21 files changed, 1427 insertions(+), 28 deletions(-) create mode 100755 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/IntArray.h create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunction.h create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyParametrization.h create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSExtrapolationState.h create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrization.h create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrizationBase.h create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSSimulationState.h create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSTruthState.h create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/IntArray.cxx create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCS1DFunction.cxx create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyParametrization.cxx create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSExtrapolationState.cxx create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrization.cxx create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrizationBase.cxx create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSSimulationState.cxx create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSTruthState.cxx diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt index f2fa6da12bc..7cc0eaa65f6 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt @@ -1,3 +1,4 @@ + ################################################################################ # Package: ISF_FastCaloSimEvent ################################################################################ @@ -11,28 +12,29 @@ atlas_depends_on_subdirs( PUBLIC Control/DataModel Control/SGTools GaudiKernel + Calorimeter/CaloGeoHelpers TileCalorimeter/TileSimEvent ) # External dependencies: find_package( CLHEP ) -find_package( ROOT COMPONENTS Core Tree MathCore Hist RIO pthread ) +find_package( ROOT COMPONENTS Core Tree MathCore Physics Hist RIO pthread ) # Remove the --as-needed linker flags: atlas_disable_as_needed() # Component(s) in the package: -atlas_add_library( ISF_FastCaloSimEvent +atlas_add_root_dictionary( ISF_FastCaloSimEvent _dictSource + ROOT_HEADERS ISF_FastCaloSimEvent/IntArray.h ) + +atlas_add_library( ISF_FastCaloSimEvent ${_dictSource} src/*.cxx PUBLIC_HEADERS ISF_FastCaloSimEvent - INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} - PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} + INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS} DEFINITIONS ${CLHEP_DEFINITIONS} - LINK_LIBRARIES ${CLHEP_LIBRARIES} DataModel SGTools GaudiKernel TileSimEvent - PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} ) + LINK_LIBRARIES ${CLHEP_LIBRARIES} ${ROOT_LIBRARIES} DataModel SGTools GaudiKernel CaloGeoHelpers TileSimEvent ) atlas_add_dictionary( ISF_FastCaloSimEventDict ISF_FastCaloSimEvent/ISF_FastCaloSimEventDict.h ISF_FastCaloSimEvent/selection.xml INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} - LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} DataModel SGTools GaudiKernel TileSimEvent ISF_FastCaloSimEvent ) - + LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} DataModel SGTools GaudiKernel CaloGeoHelpers TileSimEvent ISF_FastCaloSimEvent ) diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h new file mode 100755 index 00000000000..006e00553e4 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h @@ -0,0 +1,30 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef FastCaloSim_CaloCell_ID_h +#define FastCaloSim_CaloCell_ID_h + +#include "CaloGeoHelpers/CaloSampling.h" + +namespace CaloCell_ID_FCS { + enum CaloSample_FCS { + FirstSample=CaloSampling::PreSamplerB, + PreSamplerB=CaloSampling::PreSamplerB, EMB1=CaloSampling::EMB1, EMB2=CaloSampling::EMB2, EMB3=CaloSampling::EMB3, // LAr barrel + PreSamplerE=CaloSampling::PreSamplerE, EME1=CaloSampling::EME1, EME2=CaloSampling::EME2, EME3=CaloSampling::EME3, // LAr EM endcap + HEC0=CaloSampling::HEC0, HEC1=CaloSampling::HEC1, HEC2=CaloSampling::HEC2, HEC3=CaloSampling::HEC3, // Hadronic end cap cal. + TileBar0=CaloSampling::TileBar0, TileBar1=CaloSampling::TileBar1, TileBar2=CaloSampling::TileBar2, // Tile barrel + TileGap1=CaloSampling::TileGap1, TileGap2=CaloSampling::TileGap2, TileGap3=CaloSampling::TileGap3, // Tile gap (ITC & scint) + TileExt0=CaloSampling::TileExt0, TileExt1=CaloSampling::TileExt1, TileExt2=CaloSampling::TileExt2, // Tile extended barrel + FCAL0=CaloSampling::FCAL0, FCAL1=CaloSampling::FCAL1, FCAL2=CaloSampling::FCAL2, // Forward EM endcap + + // Beware of MiniFCAL! We don't have it, so different numbers after FCAL2 + + LastSample = CaloSampling::FCAL2, + MaxSample = LastSample+1, + noSample = -1 + }; + typedef CaloSample_FCS CaloSample; +} +#endif + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/IntArray.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/IntArray.h new file mode 100644 index 00000000000..e965dfd322a --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/IntArray.h @@ -0,0 +1,26 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef IntArray_h +#define IntArray_h + +#include "TArrayI.h" +#include "TObject.h" + +class IntArray:public TObject, public TArrayI +{ + public: + + IntArray(); + IntArray(int); + + private: + + TArrayI* m_array; + + ClassDef(IntArray,0) + +}; + +#endif diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunction.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunction.h new file mode 100644 index 00000000000..4b585de7f26 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunction.h @@ -0,0 +1,42 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TFCS1DFunction_h +#define TFCS1DFunction_h + +// STL includes +#include <string> + +#include "TH1.h" +#include "TTree.h" + +class TFCS1DFunction:public TObject +{ + public: + TFCS1DFunction(); + TFCS1DFunction(TH1* hist); + + virtual void Initialize(TH1* hist); + int testHisto(TH1* hist, std::string, float&, float&, std::string, int, std::string); + + virtual double rnd_to_fct(double rnd); + + TH1* transform(TH1*, float&, float&); + double get_maxdev(TH1* h_input, TH1* h_approx); + void tmvaregression_training(int, TTree *regTree, std::string, std::string); + double get_range_low(TH1* hist); + double tmvaregression_application(double, std::string); + static TH1* get_cumul(TH1* hist); + + private: + + //ClassDef(TFCS1DFunction,1) //TFCS1DFunction +}; + + +#if defined(__MAKECINT__) +#pragma link C++ class TFCS1DFunction; +#endif + +#endif diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyParametrization.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyParametrization.h new file mode 100644 index 00000000000..7f844be747a --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyParametrization.h @@ -0,0 +1,29 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TFCSEnergyParametrization_h +#define TFCSEnergyParametrization_h + +#include "ISF_FastCaloSimEvent/TFCSParametrization.h" + +class TFCSEnergyParametrization:public TFCSParametrization { +public: + TFCSEnergyParametrization(const char* name=0, const char* title=0); + + virtual bool is_match_Ekin_bin(int /*Ekin_bin*/) const {return true;}; + virtual bool is_match_calosample(int /*calosample*/) const {return true;}; + + // return number of energy parametrization bins + virtual int n_bins() {return 0;}; + +private: + + //ClassDef(TFCSEnergyParametrization,1) //TFCSEnergyParametrization +}; + +#if defined(__MAKECINT__) +#pragma link C++ class TFCSEnergyParametrization; +#endif + +#endif diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSExtrapolationState.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSExtrapolationState.h new file mode 100644 index 00000000000..b7d94d8abf7 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSExtrapolationState.h @@ -0,0 +1,65 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TFCSExtrapolationState_h +#define TFCSExtrapolationState_h + +#include <TObject.h> +#include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h" + +class TFCSExtrapolationState:public TObject { + public: + TFCSExtrapolationState(); + + void clear(); + + enum SUBPOS { SUBPOS_MID = 0, SUBPOS_ENT = 1, SUBPOS_EXT = 2}; //MID=middle, ENT=entrance, EXT=exit of cal layer + + void set_OK (int layer,int subpos,bool val) {m_CaloOK[layer][subpos]=val;}; + void set_eta(int layer,int subpos,double val) {m_etaCalo[layer][subpos]=val;}; + void set_phi(int layer,int subpos,double val) {m_phiCalo[layer][subpos]=val;}; + void set_r (int layer,int subpos,double val) {m_rCalo[layer][subpos]=val;}; + void set_z (int layer,int subpos,double val) {m_zCalo[layer][subpos]=val;}; + void set_d (int layer,int subpos,double val) {m_dCalo[layer][subpos]=val;}; + void set_detaBorder(int layer,int subpos,double val) {m_distetaCaloBorder[layer][subpos]=val;}; + void set_IDCaloBoundary_eta(double val) {m_IDCaloBoundary_eta=val;}; + void set_IDCaloBoundary_phi(double val) {m_IDCaloBoundary_phi=val;}; + void set_IDCaloBoundary_r(double val) {m_IDCaloBoundary_r=val;}; + void set_IDCaloBoundary_z(double val) {m_IDCaloBoundary_z=val;}; + + bool OK(int layer,int subpos) const {return m_CaloOK[layer][subpos];}; + double eta(int layer,int subpos) const {return m_etaCalo[layer][subpos];}; + double phi(int layer,int subpos) const {return m_phiCalo[layer][subpos];}; + double r(int layer,int subpos) const {return m_rCalo[layer][subpos];}; + double z(int layer,int subpos) const {return m_zCalo[layer][subpos];}; + double d(int layer,int subpos) const {return m_dCalo[layer][subpos];}; + double detaBorder(int layer,int subpos) const {return m_distetaCaloBorder[layer][subpos];}; + + double IDCaloBoundary_eta() const {return m_IDCaloBoundary_eta;}; + double IDCaloBoundary_phi() const {return m_IDCaloBoundary_phi;}; + double IDCaloBoundary_r() const {return m_IDCaloBoundary_r;}; + double IDCaloBoundary_z() const {return m_IDCaloBoundary_z;}; + + private: + bool m_CaloOK[CaloCell_ID_FCS::MaxSample][3]; + double m_etaCalo[CaloCell_ID_FCS::MaxSample][3]; + double m_phiCalo[CaloCell_ID_FCS::MaxSample][3]; + double m_rCalo[CaloCell_ID_FCS::MaxSample][3]; + double m_zCalo[CaloCell_ID_FCS::MaxSample][3]; + double m_dCalo[CaloCell_ID_FCS::MaxSample][3]; + double m_distetaCaloBorder[CaloCell_ID_FCS::MaxSample][3]; + + double m_IDCaloBoundary_eta; + double m_IDCaloBoundary_phi; + double m_IDCaloBoundary_r; + double m_IDCaloBoundary_z; + + //ClassDef(TFCSExtrapolationState,1) //TFCSExtrapolationState +}; + +#if defined(__MAKECINT__) +#pragma link C++ class TFCSExtrapolationState; +#endif + +#endif diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h new file mode 100644 index 00000000000..c8621b6cb94 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h @@ -0,0 +1,56 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TFCSPCAEnergyParametrization_h +#define TFCSPCAEnergyParametrization_h + +#include "ISF_FastCaloSimEvent/TFCSEnergyParametrization.h" +#include "ISF_FastCaloSimEvent/IntArray.h" +#include "ISF_FastCaloSimEvent/TFCS1DFunction.h" +#include "TMatrixF.h" +#include "TMatrixDSym.h" +#include "TVectorF.h" +#include "TFile.h" + +class TFCSPCAEnergyParametrization:public TFCSEnergyParametrization +{ + public: + TFCSPCAEnergyParametrization(const char* name=0, const char* title=0); + + // energies in calo layers should be returned in simulstate + virtual void simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol); + + //int n_bins() {return 0;}; + // CHANGE: return number of energy parametrization bins + int n_bins() { return (m_RelevantLayers->GetSize()+1); }; + void P2X(TMatrixD, int gNVariables, double *p, double *x, int nTest); + void P2X(int gNVariables, double *p, double *x, int nTest); + void loadInputs(TFile* file, int); + void loadInputs(TFile* file); + double InvCumulant(TH1D* hist,double y); + + private: + // PCA Matrix and NN mapping information should be stored as private member variables here + + TMatrixDSym* m_symCov; + TMatrixD* m_EV; + TVectorD* m_MeanValues; + TVectorD* m_SigmaValues; + TVectorD* m_Gauss_means; + TVectorD* m_Gauss_rms; //Gauss nur einmal speichern! + IntArray* m_RelevantLayers; + TVectorD* m_LowerBounds; + std::vector<TFCS1DFunction*> m_cumulative; + std::vector<TH1D*> h_cumulative; + + //ClassDef(TFCSPCAEnergyParametrization,1) //TFCSPCAEnergyParametrization + +}; + +#if defined(__MAKECINT__) +#pragma link C++ class TFCSPCAEnergyParametrization; +#endif + +#endif + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrization.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrization.h new file mode 100644 index 00000000000..dae0718a5a2 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrization.h @@ -0,0 +1,49 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TFCSParametrization_h +#define TFCSParametrization_h + +#include "ISF_FastCaloSimEvent/TFCSParametrizationBase.h" + +class TFCSParametrization:public ::TFCSParametrizationBase { +public: + TFCSParametrization(const char* name=0, const char* title=0); + + virtual bool is_match_pdgid(int id) const {return m_pdgid.find(id)!=m_pdgid.end();}; + virtual bool is_match_Ekin(float Ekin) const {return (Ekin>=m_Ekin_min) && (Ekin<m_Ekin_max);}; + virtual bool is_match_eta(float eta) const {return (eta>=m_eta_min) && (eta<m_eta_max);}; + + const std::set< int > &pdgid() const {return m_pdgid;}; + double Ekin_nominal() const {return m_Ekin_nominal;}; + double Ekin_min() const {return m_Ekin_min;}; + double Ekin_max() const {return m_Ekin_max;}; + double eta_nominal() const {return m_eta_nominal;}; + double eta_min() const {return m_eta_min;}; + double eta_max() const {return m_eta_max;}; + + void set_pdgid(int id); + void add_pdgid(int id); + void clear_pdgid(); + + void set_Ekin_nominal(double min); + void set_Ekin_min(double min); + void set_Ekin_max(double max); + void set_eta_nominal(double min); + void set_eta_min(double min); + void set_eta_max(double max); + +private: + std::set< int > m_pdgid; + double m_Ekin_nominal,m_Ekin_min,m_Ekin_max; + double m_eta_nominal,m_eta_min,m_eta_max; + + //ClassDef(TFCSParametrization,1) //TFCSParametrization +}; + +#if defined(__MAKECINT__) +#pragma link C++ class TFCSParametrization; +#endif + +#endif diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrizationBase.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrizationBase.h new file mode 100644 index 00000000000..22bb448c6d1 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrizationBase.h @@ -0,0 +1,49 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TFCSParametrizationBase_h +#define TFCSParametrizationBase_h + +#include <TNamed.h> +#include <set> +#include "ISF_FastCaloSimEvent/TFCSSimulationState.h" +#include "ISF_FastCaloSimEvent/TFCSTruthState.h" +#include "ISF_FastCaloSimEvent/TFCSExtrapolationState.h" + +class TFCSParametrizationBase:public TNamed { +public: + TFCSParametrizationBase(const char* name=0, const char* title=0); + + virtual bool is_match_pdgid(int /*id*/) const {return false;}; + virtual bool is_match_Ekin(float /*Ekin*/) const {return false;}; + virtual bool is_match_eta(float /*eta*/) const {return false;}; + + virtual bool is_match_Ekin_bin(int /*Ekin_bin*/) const {return false;}; + virtual bool is_match_calosample(int /*calosample*/) const {return false;}; + + virtual const std::set< int > &pdgid() const {return m_no_pdgid;}; + virtual double Ekin_nominal() const {return 0;}; + virtual double Ekin_min() const {return 0;}; + virtual double Ekin_max() const {return 0;}; + virtual double eta_nominal() const {return 100;}; + virtual double eta_min() const {return 100;}; + virtual double eta_max() const {return 100;}; + + // Do some simulation. Result should be returned in simulstate + // Simulate all energies in calo layers for energy parametrizations + // Simulate one HIT for later shape parametrizations (TO BE DISCUSSED!) + virtual void simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol); + + void Print(Option_t *option = "") const; +private: + static std::set< int > m_no_pdgid; + + //ClassDef(TFCSParametrizationBase,1) //TFCSParametrizationBase +}; + +#if defined(__MAKECINT__) +#pragma link C++ class TFCSParametrizationBase; +#endif + +#endif diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSSimulationState.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSSimulationState.h new file mode 100644 index 00000000000..dbaeaf7fb9f --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSSimulationState.h @@ -0,0 +1,43 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TFCSSimulationState_h +#define TFCSSimulationState_h + +#include <TObject.h> +#include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h" + +class TFCSSimulationState:public TObject +{ + public: + TFCSSimulationState(); + + bool is_valid() const {return m_Ebin>=0;}; + double E() const {return m_Etot;}; + double E(int sample) const {return m_E[sample];}; + int Ebin() const {return m_Ebin;}; + + void set_Ebin(int bin) {m_Ebin=bin;}; + void set_E(int sample,double Esample) { m_E[sample]=Esample; } ; + void set_E(double E) { m_Etot=E; } ; + void add_E(int sample,double Esample) { m_E[sample]+=Esample;m_Etot+=Esample; }; + + //empty function so far + //not sure if we should keep the function here or rather write a deposit_cell function or similar + virtual void deposit_HIT(int sample,double hit_eta,double hit_phi,double hit_weight); + + void clear(); + private: + int m_Ebin; + double m_Etot; + double m_E[CaloCell_ID_FCS::MaxSample]; + + //ClassDef(TFCSSimulationState,1) //TFCSSimulationState +}; + +#if defined(__MAKECINT__) +#pragma link C++ class TFCSSimulationState; +#endif + +#endif diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSTruthState.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSTruthState.h new file mode 100644 index 00000000000..b56e6a0aada --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSTruthState.h @@ -0,0 +1,33 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TFCSTruthState_h +#define TFCSTruthState_h + +#include <TLorentzVector.h> + +class TFCSTruthState:public TLorentzVector { + public: + TFCSTruthState(); + TFCSTruthState(Double_t x, Double_t y, Double_t z, Double_t t, int pdgid); + + void set_pdgid(int val) {m_pdgid=val;}; + void set_vertex(const TLorentzVector& val) {m_vertex=val;}; + void set_vertex(Double_t x, Double_t y, Double_t z, Double_t t=0) {m_vertex.SetXYZT(x,y,z,t);}; + + int pdgid() const {return m_pdgid;}; + double Ekin() const {return E()-M();}; + const TLorentzVector& vertex() const {return m_vertex;}; + private: + int m_pdgid; + TLorentzVector m_vertex; + + //ClassDef(TFCSTruthState,1) //TFCSTruthState +}; + +#if defined(__MAKECINT__) +#pragma link C++ class TFCSTruthState; +#endif + +#endif diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/cmt/requirements b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/cmt/requirements index 1829e7bfc6f..2c2a60f97d2 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/cmt/requirements +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/cmt/requirements @@ -1,27 +1,19 @@ package ISF_FastCaloSimEvent public -use AtlasPolicy AtlasPolicy-* -use AtlasCLHEP AtlasCLHEP-* External -use GaudiInterface GaudiInterface-* External -#use LArSimEvent LArSimEvent-* LArCalorimeter/ -use TileSimEvent TileSimEvent-* TileCalorimeter -use CLIDSvc CLIDSvc-* Control -use DataModel DataModel-* Control +use AtlasPolicy AtlasPolicy-* +use AtlasCLHEP AtlasCLHEP-* External +use GaudiInterface GaudiInterface-* External +use TileSimEvent TileSimEvent-* TileCalorimeter +use CLIDSvc CLIDSvc-* Control +use DataModel DataModel-* Control +use AtlasROOT AtlasROOT-* External +use CaloGeoHelpers CaloGeoHelpers-* Calorimeter -#use AthenaBaseComps AthenaBaseComps-* Control -#use AtlasHepMC AtlasHepMC-* External -#use BarcodeInterfaces BarcodeInterfaces-* Simulation/Barcode -#use AtlasDetDescr AtlasDetDescr-* DetectorDescription -#use GeoPrimitives GeoPrimitives-* DetectorDescription -end_public - -private - -public - -library ISF_FastCaloSimEvent *.cxx -apply_pattern installed_library +library ISF_FastCaloSimEventLib *.cxx +#apply_pattern installed_library +apply_pattern named_installed_library library=ISF_FastCaloSimEventLib +apply_pattern have_root_headers root_headers="IntArray.h" headers_lib="ISF_FastCaloSimEventLib" #for dictionaries private diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/IntArray.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/IntArray.cxx new file mode 100644 index 00000000000..43a1ea7e978 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/IntArray.cxx @@ -0,0 +1,31 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "ISF_FastCaloSimEvent/IntArray.h" + +#include "TArrayI.h" + +//========================== +//======= IntArray ========= +//========================== + +using namespace std; + +IntArray::IntArray(int n) +{ + m_array=new TArrayI(n); +} + +IntArray::IntArray() +{ + m_array=new TArrayI(); +} + + +//============================================= +//========== ROOT persistency stuff =========== +//============================================= + +ClassImp(IntArray) + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCS1DFunction.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCS1DFunction.cxx new file mode 100644 index 00000000000..fa4df6975c5 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCS1DFunction.cxx @@ -0,0 +1,371 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "ISF_FastCaloSimEvent/TFCS1DFunction.h" + +#include "TMVA/Config.h" +#include "TMVA/Tools.h" +#include "TMVA/Reader.h" +#include "TMVA/Factory.h" +#include "TFile.h" + +#include "TRandom1.h" + + +//============================================= +//======= TFCS1DFunction ========= +//============================================= + +TFCS1DFunction::TFCS1DFunction() +{ +} + +TFCS1DFunction::TFCS1DFunction(TH1* hist) +{ + Initialize(hist); +} + +void TFCS1DFunction::Initialize(TH1* /*hist*/) +{ +} + +int TFCS1DFunction::testHisto(TH1* hist, std::string weightfilename, float &rangeval, float &startval, std::string outfilename, int skip_regression, std::string label) +{ + + //transform the histogram + TH1* h_transf=transform(hist, rangeval, startval); h_transf->SetName("h_transf"); + + //Turn the histogram into a tree: + std::vector<double> contents; + std::vector<double> centers; + for(int b=1;b<=h_transf->GetNbinsX();b++) + { + contents.push_back(h_transf->GetBinContent(b)); + centers.push_back(h_transf->GetBinCenter(b)); + } + + TTree* tree=new TTree("tree","tree"); + Float_t x,y; + tree->Branch("x",&x,"x/F"); + tree->Branch("y",&y,"y/F"); + + for(unsigned int i=0;i<centers.size();i++) + { + y=(Float_t)(contents[i]); //xvals are the BinContents + x=(Float_t)(centers[i]); //yvals are the BinCenters + + tree->Fill(); + } + + double range_low=get_range_low(h_transf); + + TRandom1* myRandom=new TRandom1(); myRandom->SetSeed(0); + int Ntoys=10000; + + int do_iteration=1; + int do_range=1; + double maxdev=100; + int neurons=2; + + double maxdev_cut=5; //5!! + + std::vector<TH1*> histos; + while(!skip_regression && maxdev>maxdev_cut && neurons<11) //11 + { + + tmvaregression_training(neurons, tree, weightfilename, outfilename); + + std::cout<<"Testing the regression with "<<Ntoys<<" toys"<<std::endl; + TH1* h_output=(TH1*)h_transf->Clone("h_output"); + h_output->Reset(); + for(int i=0;i<Ntoys;i++) + { + double random=myRandom->Uniform(1); + if(do_range && random<range_low) random=range_low; + double value =tmvaregression_application(random,weightfilename); + h_output->Fill(value); + } + + TH1* h_cumul=get_cumul(h_output); + h_cumul->SetName(Form("h_cumul_neurons%i",neurons)); + histos.push_back(h_cumul); + + maxdev=get_maxdev(h_transf,h_cumul); + std::cout<<"---> Neurons="<<neurons<<" MAXDEV="<<maxdev<<"%"<<std::endl; + neurons++; + if(!do_iteration) break; + } + + TH1* histclone=(TH1*)hist->Clone("histclone"); + + TFile* out_iteration=new TFile(Form("output/iteration_%s.root",label.c_str()),"RECREATE"); + for(unsigned h=0;h<histos.size();h++) + { + out_iteration->Add(histos[h]); + } + out_iteration->Add(histclone); + out_iteration->Write(); + out_iteration->Close(); + + int regression_success=1; + if(maxdev>maxdev_cut) regression_success=0; + + int status=0; + if(regression_success) + { + std::cout<<"Regression successful. Weights are stored."<<std::endl; + if(rangeval<0) status=1; + if(rangeval>=0) status=2; + } + + + if(!regression_success) + { + std::cout<<"Regression failed. Histogram is stored."<<std::endl; + status=3; + } //!success + + return status; + +} + +double TFCS1DFunction::rnd_to_fct(double /*rnd*/) +{ + + return 0; +} + + +double TFCS1DFunction::get_range_low(TH1* hist) +{ + double range_low=0.0; + int bin_start=-1; + for(int b=1;b<=hist->GetNbinsX();b++) + { + if(hist->GetBinContent(b)>0 && bin_start<0) + { + bin_start=b; + range_low=hist->GetBinContent(b); + b=hist->GetNbinsX()+1; + } + } + return range_low; +} + +void TFCS1DFunction::tmvaregression_training(int neurons, TTree *regTree, std::string weightfile, std::string outfilename) +{ + using namespace TMVA; + + TString myMethodList = "" ; + TMVA::Tools::Instance(); + std::map<std::string,int> Use; + + Use["PDERS"] = 0; Use["PDEFoam"] = 0; Use["KNN"] = 0; Use["LD"] = 0; Use["FDA_GA"] = 0; Use["FDA_MC"] = 0; + Use["FDA_MT"] = 0; Use["FDA_GAMT"] = 0; Use["MLP"] = 1; Use["SVM"] = 0; Use["BDT"] = 0; Use["BDTG"] = 0; + + std::cout << std::endl; std::cout << "==> Start TMVARegression with "<<neurons<<" Neurons "<<std::endl; + + if(myMethodList != "") + { + for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) it->second = 0; + std::vector<TString> mlist = gTools().SplitString( myMethodList, ',' ); + for (UInt_t i=0; i<mlist.size(); i++) + { + std::string regMethod(mlist[i]); + if (Use.find(regMethod) == Use.end()) + { + std::cout << "Method \"" << regMethod << "\" not known in TMVA under this name. Choose among the following:" << std::endl; + for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) std::cout << it->first << " "; + std::cout << std::endl; + return; + } + Use[regMethod] = 1; + } + } + +// TString outfileName("TMVAReg.root"); + TFile* outputFile = TFile::Open( outfilename.c_str(), "RECREATE" ); + TMVA::Factory *factory = new TMVA::Factory( "TMVARegression", outputFile, "!V:!Silent:Color:DrawProgressBar" ); + TString dirname=Form("%s/",weightfile.c_str()); + (TMVA::gConfig().GetIONames()).fWeightFileDir = dirname; + factory->AddVariable( "y", "y", 'F' ); + factory->AddTarget( "x" ); + //TFile *input(0); + Double_t regWeight = 1.0; + + factory->AddRegressionTree( regTree, regWeight ); + TCut mycut = ""; + factory->PrepareTrainingAndTestTree( mycut,"nTrain_Regression=0:nTest_Regression=0:SplitMode=Random:NormMode=NumEvents:!V" ); + + if(Use["MLP"]) + //factory->BookMethod( TMVA::Types::kMLP, "MLP", Form("!H:!V:VarTransform=Norm:NeuronType=tanh:NCycles=20000:HiddenLayers=%i:TestRate=6:TrainingMethod=BFGS:Sampling=0.3:SamplingEpoch=0.8:ConvergenceImprove=1e-6:ConvergenceTests=15:!UseRegulator",neurons) ); + factory->BookMethod( TMVA::Types::kMLP, "MLP", Form("!H:!V:NeuronType=tanh:NCycles=20000:HiddenLayers=%i:TestRate=6:TrainingMethod=BFGS:Sampling=0.3:SamplingEpoch=0.8:ConvergenceImprove=1e-6:ConvergenceTests=15:!UseRegulator",neurons) ); + + // Train MVAs using the set of training events + factory->TrainAllMethods(); + + // ---- Evaluate all MVAs using the set of test events + factory->TestAllMethods(); + + // ----- Evaluate and compare performance of all configured MVAs + factory->EvaluateAllMethods(); + + // Save the output + outputFile->Close(); + + std::cout << "==> Wrote root file: " << outputFile->GetName() << std::endl; + std::cout << "==> TMVARegression is done!" << std::endl; + + delete factory; + +} + + +TH1* TFCS1DFunction::transform(TH1* h_input, float &rangeval, float& startval) +{ + + bool do_transform=false; + float xmin=h_input->GetXaxis()->GetXmin(); + float xmax=h_input->GetXaxis()->GetXmax(); + if(xmin<0 || xmax>1) do_transform=true; + + TH1D* h_out; + + if(do_transform) + { + int nbins=h_input->GetNbinsX(); + double min=0; + double max=1; + h_out=new TH1D("h_out","h_out",nbins,min,max); + + for(int b=1;b<=nbins;b++) + h_out->SetBinContent(b,h_input->GetBinContent(b)); + + //store the inital range + rangeval=xmax-xmin; + startval=xmin; + } + if(!do_transform) + { + rangeval=-1; + h_out=(TH1D*)h_input->Clone("h_out"); + } + return h_out; + +} + + +double TFCS1DFunction::get_maxdev(TH1* h_input1, TH1* h_approx1) +{ + + TH1D* h_input =(TH1D*)h_input1->Clone("h_input"); + TH1D* h_approx=(TH1D*)h_approx1->Clone("h_approx"); + + double maxdev=0.0; + + //normalize the histos to the same area: + double integral_input=h_input->Integral(); + double integral_approx=0.0; + for(int b=1;b<=h_input->GetNbinsX();b++) + integral_approx+=h_approx->GetBinContent(h_approx->FindBin(h_input->GetBinCenter(b))); + h_approx->Scale(integral_input/integral_approx); + + double ymax=h_approx->GetBinContent(h_approx->GetNbinsX())-h_approx->GetBinContent(h_approx->GetMinimumBin()); + for(int i=1;i<=h_input->GetNbinsX();i++) + { + double val=fabs(h_approx->GetBinContent(h_approx->FindBin(h_input->GetBinCenter(i)))-h_input->GetBinContent(i))/ymax; + if(val>maxdev) maxdev=val; + } + + delete h_input; + delete h_approx; + + return maxdev*100.0; + +} + + +double TFCS1DFunction::tmvaregression_application(double uniform, std::string weightfile) +{ + + using namespace TMVA; + + TString myMethodList = "" ; + TMVA::Tools::Instance(); + + std::map<std::string,int> Use; + + // --- Mutidimensional likelihood and Nearest-Neighbour methods + Use["PDERS"] = 0; Use["PDEFoam"] = 0; Use["KNN"] = 0; + Use["LD"] = 0; Use["FDA_GA"] = 0; Use["FDA_MC"] = 0; + Use["FDA_MT"] = 0; Use["FDA_GAMT"] = 0; Use["MLP"] = 1; + Use["SVM"] = 0; Use["BDT"] = 0; Use["BDTG"] = 0; + // --------------------------------------------------------------- + + // Select methods (don't look at this code - not of interest) + if(myMethodList != "") + { + for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) it->second = 0; + std::vector<TString> mlist = gTools().SplitString( myMethodList, ',' ); + for (UInt_t i=0; i<mlist.size(); i++) + { + std::string regMethod(mlist[i]); + if (Use.find(regMethod) == Use.end()) + { + std::cout << "Method \"" << regMethod << "\" not known in TMVA under this name. Choose among the following:" << std::endl; + for(std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) std::cout << it->first << " "; + std::cout << std::endl; + return 0; + } + Use[regMethod] = 1; + } + } + + // -------------------------------------------------------------------------------------------------- + + TMVA::Reader *reader = new TMVA::Reader( "!Color:Silent"); + + Float_t y=uniform; + reader->AddVariable( "y", &y ); + + TString dir = Form("%s/",weightfile.c_str()); + TString prefix = "TMVARegression"; + + // Book method(s) + for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) + { + if (it->second) + { + TString methodName = it->first + " method"; + TString weightfile = dir + prefix + "_" + TString(it->first) + ".weights.xml"; + reader->BookMVA( methodName, weightfile ); + } + } + + Float_t val = (reader->EvaluateRegression("MLP method"))[0]; + delete reader; + return val; + +} + + +TH1* TFCS1DFunction::get_cumul(TH1* hist) +{ + TH1D* h_cumul=(TH1D*)hist->Clone("h_cumul"); + double sum=0; + for(int b=1;b<=h_cumul->GetNbinsX();b++) + { + sum+=hist->GetBinContent(b); + h_cumul->SetBinContent(b,sum); + } + return h_cumul; +} + + +//============================================= +//========== ROOT persistency stuff =========== +//============================================= + +//ClassImp(TFCS1DFunction) + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyParametrization.cxx new file mode 100644 index 00000000000..8121cc70255 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyParametrization.cxx @@ -0,0 +1,21 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "ISF_FastCaloSimEvent/TFCSEnergyParametrization.h" +#include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h" + +//============================================= +//======= TFCSEnergyParametrization ========= +//============================================= + +TFCSEnergyParametrization::TFCSEnergyParametrization(const char* name, const char* title):TFCSParametrization(name,title) +{ +} + +//============================================= +//========== ROOT persistency stuff =========== +//============================================= + +//ClassImp(TFCSEnergyParametrization) + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSExtrapolationState.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSExtrapolationState.cxx new file mode 100644 index 00000000000..191f338b2ef --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSExtrapolationState.cxx @@ -0,0 +1,42 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "ISF_FastCaloSimEvent/TFCSExtrapolationState.h" + +//============================================= +//======= TFCSExtrapolationState ========= +//============================================= + +TFCSExtrapolationState::TFCSExtrapolationState() +{ + clear(); +} + +void TFCSExtrapolationState::clear() +{ + for(int i=0;i<CaloCell_ID_FCS::MaxSample;++i) { + for(int j=0;j<3;++j) { + m_CaloOK[i][j]=false; + m_etaCalo[i][j]=-999; + m_phiCalo[i][j]=-999; + m_rCalo[i][j]=0; + m_zCalo[i][j]=0; + m_dCalo[i][j]=0; + m_distetaCaloBorder[i][j]=0; + } + } + + m_IDCaloBoundary_eta=-999; + m_IDCaloBoundary_phi=-999; + m_IDCaloBoundary_r=0; + m_IDCaloBoundary_z=0; +} + + +//============================================= +//========== ROOT persistency stuff =========== +//============================================= + +//ClassImp(TFCSExtrapolationState) + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx new file mode 100644 index 00000000000..d198e6cd36b --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx @@ -0,0 +1,352 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h" +#include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h" + +#include "TFile.h" +#include <iostream> +#include "TKey.h" +#include "TClass.h" +#include "TRandom3.h" +#include "TMatrixD.h" +#include "TMatrixDSymEigen.h" + +//============================================= +//======= TFCSPCAEnergyParametrization ========= +//============================================= + +TFCSPCAEnergyParametrization::TFCSPCAEnergyParametrization(const char* name, const char* title):TFCSEnergyParametrization(name,title) +{ +} + +void TFCSPCAEnergyParametrization::simulate(TFCSSimulationState& simulstate,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* /*extrapol*/) +{ + + int verbose=0; + + int thisbin=simulstate.Ebin(); + if(verbose) std::cout<<"--- Simulation of bin "<<thisbin<<std::endl; + + /* + for(unsigned int i=0;i<CaloCell_ID_FCS::MaxSample;++i) + simulstate.add_E(i,truth->Ekin()/CaloCell_ID_FCS::MaxSample); + */ + + TRandom3* Random3=new TRandom3(); + Random3->SetSeed(0); + + //get relevant layers for that bin + std::vector<std::string> layer; + std::vector<int> layerNr; + + for(int i=0;i<m_RelevantLayers->GetSize();i++) + layerNr.push_back(m_RelevantLayers->GetAt(i)); + + for(unsigned int i=0;i<layerNr.size();i++) + { + std::string thislayer=Form("layer%i",layerNr[i]); + layer.push_back(thislayer); + } + layer.push_back("totalE"); + + //for(unsigned int i=0;i<layer.size();i++) cout<<layer[i]<<endl; + + double* gauss_means=m_Gauss_means->GetMatrixArray(); + double* gauss_rms =m_Gauss_rms->GetMatrixArray(); + double* lowerBounds=m_LowerBounds->GetMatrixArray(); + + double output_data[layer.size()]; + double input_data[layer.size()]; + + for(unsigned int l=0;l<layer.size();l++) + { + double mean=gauss_means[l]; + double rms =gauss_rms[l]; + double gauszz=Random3->Gaus(mean,rms); + input_data[l]=gauszz; + } + + TMatrixDSymEigen* eigen=new TMatrixDSymEigen(*m_symCov); + TMatrixD myEvectors(eigen->GetEigenVectors()); + + //P2X(myEvectors, layer.size(), input_data, output_data, layer.size()); + P2X(layer.size(), input_data, output_data, layer.size()); + + double simdata_uniform[layer.size()]; + double simdata[layer.size()]; + double simdata_scaled[layer.size()]; + double sum_fraction=0.0; //this is needed for the rescaling in the next step + for(unsigned int l=0;l<layer.size();l++) + { + simdata_uniform[l]=(TMath::Erf(output_data[l]/1.414213562)+1)/2.f; + + if(simdata_uniform[l]<lowerBounds[l]) simdata_uniform[l]=lowerBounds[l]; + + simdata[l]=m_cumulative[l]->rnd_to_fct(simdata_uniform[l]); + //simdata[l]=InvCumulant(h_cumulative[l],simdata_uniform[l]); + + if(simdata[l]<0) simdata[l]=0; + if(layer[l]!="totalE" && simdata[l]>1) simdata[l]=1; + if(layer[l]!="totalE") sum_fraction+=simdata[l]; + } + + //RECSCALE SO THAT THE SUM GIVES 1 + double scale=1.0/sum_fraction; + + sum_fraction=0.0; + + for(unsigned int l=0;l<layer.size();l++) + { + if(layer[l]!="totalE") + { + simdata_scaled[l]=simdata[l]*scale; + if(l<layerNr.size()) + { + sum_fraction+=simdata_scaled[l]; + } + } + } //for layer + + //save the result in the SimulState object: + + for(int s=0;s<CaloCell_ID_FCS::MaxSample;s++) + { + double energy=0.0; + for(unsigned int l=0;l<layerNr.size();l++) + { + if(layerNr[l]==s) + energy=simdata_scaled[l]; + } + simulstate.set_E(s,energy); + } + simulstate.set_E(simdata[layer.size()-1]); + if(verbose) + { + std::cout<<"Result:"<<std::endl; + for(unsigned int l=0;l<layerNr.size();l++) + std::cout<<" Energy fraction deposited in Layer "<<layerNr[l]<<" = "<<simdata_scaled[l]<<std::endl; + std::cout<<" Sum of energy fractions = "<<sum_fraction<<std::endl; + std::cout<<" Total Energy = "<<simulstate.E()<<" MeV "<<std::endl; + } + + delete Random3; + +} + + +void TFCSPCAEnergyParametrization::P2X(TMatrixD EigenVectors, int gNVariables, double *p, double *x, int nTest) +{ + + double* gSigmaValues = m_SigmaValues->GetMatrixArray(); + double* gMeanValues = m_MeanValues->GetMatrixArray(); + double* gEigenVectors = EigenVectors.GetMatrixArray(); + + for(int i = 0; i < gNVariables; i++) + { + x[i] = gMeanValues[i]; + for(int j = 0; j < nTest; j++) + { + x[i] += p[j] * gSigmaValues[i] * (double)(gEigenVectors[i * gNVariables + j]); + } + } +} + +void TFCSPCAEnergyParametrization::P2X(int gNVariables, double *p, double *x, int nTest) +{ + + double* gSigmaValues = m_SigmaValues->GetMatrixArray(); + double* gMeanValues = m_MeanValues->GetMatrixArray(); + double* gEigenVectors = m_EV->GetMatrixArray(); + + for(int i = 0; i < gNVariables; i++) + { + x[i] = gMeanValues[i]; + for(int j = 0; j < nTest; j++) + { + x[i] += p[j] * gSigmaValues[i] * (double)(gEigenVectors[i * gNVariables + j]); + } + } +} + +double TFCSPCAEnergyParametrization::InvCumulant(TH1D* hist,double y) +{ + int bin = 0; + int nbin = hist->GetNbinsX(); + double min = 99999999; + for (int i=1; i<=hist->GetNbinsX()-2; i++) + { + if(fabs(hist->GetBinContent(i)-y)<min) + { + min = fabs(hist->GetBinContent(i)-y); + bin = i ; + } + } + bin = TMath::Max(bin,1); + bin = TMath::Min(bin,hist->GetNbinsX()); + + //std::cout<<bin <<std::endl; + + double AvNuEvPerBin; + double Tampon = 0 ; + for (int i=1; i<=nbin; i++) + { + Tampon += hist->GetBinContent(i); + } + + AvNuEvPerBin = Tampon/nbin; + + double x; + double x0, x1, y0, y1; + double total = hist->GetNbinsX()*AvNuEvPerBin; + double supmin = 0.5/total; + + x0 = hist->GetBinLowEdge(TMath::Max(bin,1)); + x1 = hist->GetBinLowEdge(TMath::Min(bin,hist->GetNbinsX())+1); + + y0 = hist->GetBinContent(TMath::Max(bin-1,0)); // Y0 = F(x0); Y0 >= 0 + y1 = hist->GetBinContent(TMath::Min(bin, hist->GetNbinsX()+1)); // Y1 = F(x1); Y1 <= 1 + + //Zero bin + if (bin == 0) + { + y0 = supmin; + y1 = supmin; + } + if (bin == 1) { + y0 = supmin; + } + if (bin > hist->GetNbinsX()) { + y0 = 1.-supmin; + y1 = 1.-supmin; + } + if (bin == hist->GetNbinsX()) { + y1 = 1.-supmin; + } + + //////////////////////// + + if(y0 == x1) + { + x = x0; + } + else + { + x = x0 + (y-y0)*(x1-x0)/(y1-y0); + } + + return x; + +} + + +void TFCSPCAEnergyParametrization::loadInputs(TFile* file, int bin) +{ + + file->cd(Form("bin%i",bin)); + + m_symCov =(TMatrixDSym*)gDirectory->Get("symCov"); + m_EV =(TMatrixD*)gDirectory->Get("EigenVectors"); + m_MeanValues =(TVectorD*)gDirectory->Get("MeanValues"); + m_SigmaValues =(TVectorD*)gDirectory->Get("SigmaValues"); + m_Gauss_means =(TVectorD*)gDirectory->Get("Gauss_means"); + m_Gauss_rms =(TVectorD*)gDirectory->Get("Gauss_rms"); + m_LowerBounds =(TVectorD*)gDirectory->Get("LowerBounds"); + m_RelevantLayers =(IntArray*)gDirectory->Get("RelevantLayers"); + + std::vector<std::string> layer; + std::vector<int> layerNr; + for(int i=0;i<m_RelevantLayers->GetSize();i++) + layerNr.push_back(m_RelevantLayers->GetAt(i)); + for(unsigned int i=0;i<layerNr.size();i++) + { + std::string thislayer=Form("layer%i",layerNr[i]); + layer.push_back(thislayer); + } + layer.push_back("totalE"); + + //for(unsigned int i=0;i<layer.size();i++) cout<<layer[i]<<endl; + + for(unsigned int l=0;l<layer.size();l++) + { + file->cd(Form("bin%i/%s",bin,layer[l].c_str())); + TFCS1DFunction* fct; + fct=(TFCS1DFunction*)gDirectory->Get("TFCS1DFunctionRegression"); + if(!fct) + { + fct=(TFCS1DFunction*)gDirectory->Get("TFCS1DFunctionRegressionTF"); + } + if(!fct) + { + fct=(TFCS1DFunction*)gDirectory->Get("TFCS1DFunctionHistogram"); + } + m_cumulative.push_back(fct); + + file->cd(Form("bin%i",bin)); + TH1D* hist=(TH1D*)gDirectory->Get(Form("h_cumul_%s",layer[l].c_str())); hist->SetName(Form("h_cumul_%s",layer[l].c_str())); + h_cumulative.push_back(hist); + + } + +} + +void TFCSPCAEnergyParametrization::loadInputs(TFile* file) +{ + + file->cd("pca"); + + m_symCov =(TMatrixDSym*)gDirectory->Get("symCov"); + m_EV =(TMatrixD*)gDirectory->Get("EigenVectors"); + m_MeanValues =(TVectorD*)gDirectory->Get("MeanValues"); + m_SigmaValues =(TVectorD*)gDirectory->Get("SigmaValues"); + m_Gauss_means =(TVectorD*)gDirectory->Get("Gauss_means"); + m_Gauss_rms =(TVectorD*)gDirectory->Get("Gauss_rms"); + m_LowerBounds =(TVectorD*)gDirectory->Get("LowerBounds"); + m_RelevantLayers =(IntArray*)gDirectory->Get("RelevantLayers"); + + std::vector<std::string> layer; + std::vector<int> layerNr; + for(int i=0;i<m_RelevantLayers->GetSize();i++) + layerNr.push_back(m_RelevantLayers->GetAt(i)); + for(unsigned int i=0;i<layerNr.size();i++) + { + std::string thislayer=Form("layer%i",layerNr[i]); + layer.push_back(thislayer); + } + layer.push_back("totalE"); + + //for(unsigned int i=0;i<layer.size();i++) cout<<layer[i]<<endl; + + for(unsigned int l=0;l<layer.size();l++) + { + file->cd(Form("%s",layer[l].c_str())); + TFCS1DFunction* fct; + fct=(TFCS1DFunction*)gDirectory->Get("TFCS1DFunctionRegression"); + if(!fct) + { + fct=(TFCS1DFunction*)gDirectory->Get("TFCS1DFunctionRegressionTF"); + } + if(!fct) + { + fct=(TFCS1DFunction*)gDirectory->Get("TFCS1DFunctionHistogram"); + } + m_cumulative.push_back(fct); + + } + + file->cd("pca"); + for(unsigned int l=0;l<layer.size();l++) + { + TH1D* hist=(TH1D*)gDirectory->Get(Form("h_cumul_%s",layer[l].c_str())); hist->SetName(Form("h_cumul_%s",layer[l].c_str())); + h_cumulative.push_back(hist); + } + +} + +//============================================= +//========== ROOT persistency stuff =========== +//============================================= + +//ClassImp(TFCSPCAEnergyParametrization) + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrization.cxx new file mode 100644 index 00000000000..1b3453aa9cd --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrization.cxx @@ -0,0 +1,67 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "ISF_FastCaloSimEvent/TFCSParametrization.h" + +//============================================= +//======= TFCSParametrization ========= +//============================================= + +TFCSParametrization::TFCSParametrization(const char* name, const char* title):TFCSParametrizationBase(name,title),m_Ekin_nominal(0),m_Ekin_min(0),m_Ekin_max(0),m_eta_nominal(100),m_eta_min(100),m_eta_max(100) +{ +} + +void TFCSParametrization::set_pdgid(int id) +{ + m_pdgid.clear(); + m_pdgid.insert(id); +} + +void TFCSParametrization::add_pdgid(int id) +{ + m_pdgid.insert(id); +} + +void TFCSParametrization::clear_pdgid() +{ + m_pdgid.clear(); +} + +void TFCSParametrization::set_Ekin_nominal(double nominal) +{ + m_Ekin_nominal=nominal; +} + +void TFCSParametrization::set_Ekin_min(double min) +{ + m_Ekin_min=min; +} + +void TFCSParametrization::set_Ekin_max(double max) +{ + m_Ekin_max=max; +} + + +void TFCSParametrization::set_eta_nominal(double nominal) +{ + m_eta_nominal=nominal; +} + +void TFCSParametrization::set_eta_min(double min) +{ + m_eta_min=min; +} + +void TFCSParametrization::set_eta_max(double max) +{ + m_eta_max=max; +} + +//============================================= +//========== ROOT persistency stuff =========== +//============================================= + +//ClassImp(TFCSParametrization) + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrizationBase.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrizationBase.cxx new file mode 100644 index 00000000000..28a3180a3da --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrizationBase.cxx @@ -0,0 +1,39 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "ISF_FastCaloSimEvent/TFCSParametrizationBase.h" +#include <iostream> + +//============================================= +//======= TFCSParametrizationBase ========= +//============================================= + +std::set< int > TFCSParametrizationBase::m_no_pdgid; + +TFCSParametrizationBase::TFCSParametrizationBase(const char* name, const char* title):TNamed(name,title) +{ +} + +void TFCSParametrizationBase::simulate(TFCSSimulationState& /*simulstate*/,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* /*extrapol*/) +{ +} + +void TFCSParametrizationBase::Print(Option_t *option) const +{ + TNamed::Print(option); + + std::cout <<" PDGID: "; + for (std::set<int>::iterator it=pdgid().begin(); it!=pdgid().end(); ++it) std::cout << *it << ", "; + std::cout << std::endl; + + std::cout <<" Ekin="<<Ekin_nominal()<<" MeV, range "<<Ekin_min()<<" MeV < Ekin < "<<Ekin_max()<<" MeV"<<std::endl; + std::cout <<" eta="<<eta_nominal()<<", range "<<eta_min()<<" < eta < "<<eta_max()<<std::endl; +} + +//============================================= +//========== ROOT persistency stuff =========== +//============================================= + +//ClassImp(TFCSParametrizationBase) + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSSimulationState.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSSimulationState.cxx new file mode 100644 index 00000000000..c5669b9f609 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSSimulationState.cxx @@ -0,0 +1,35 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "ISF_FastCaloSimEvent/TFCSSimulationState.h" + +//============================================= +//======= TFCSSimulationState ========= +//============================================= + +TFCSSimulationState::TFCSSimulationState() +{ + clear(); +} + +void TFCSSimulationState::clear() +{ + m_Ebin=-1; + m_Etot=0; + for(int i=0;i<CaloCell_ID_FCS::MaxSample;++i) + { + m_E[i]=0; + } +} + +void TFCSSimulationState::deposit_HIT(int /*sample*/,double /*hit_eta*/,double /*hit_phi*/,double /*hit_weight*/) +{ +} + +//============================================= +//========== ROOT persistency stuff =========== +//============================================= + +//ClassImp(TFCSSimulationState) + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSTruthState.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSTruthState.cxx new file mode 100644 index 00000000000..272556ad143 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSTruthState.cxx @@ -0,0 +1,25 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "ISF_FastCaloSimEvent/TFCSTruthState.h" + +//============================================= +//======= TFCSTruthState ========= +//============================================= + +TFCSTruthState::TFCSTruthState():TLorentzVector(),m_pdgid(0),m_vertex(0,0,0,0) +{ +} + +TFCSTruthState::TFCSTruthState(Double_t x, Double_t y, Double_t z, Double_t t, int pdgid):TLorentzVector(x,y,z,t),m_vertex(0,0,0,0) +{ + m_pdgid=pdgid; +} + +//============================================= +//========== ROOT persistency stuff =========== +//============================================= + +//ClassImp(TFCSTruthState) + -- GitLab