From dd0cab83345e00406e33ce4604c516f5cf8fcc45 Mon Sep 17 00:00:00 2001 From: Joshua Falco Beirer Date: Tue, 10 Nov 2020 16:44:17 +0100 Subject: [PATCH 1/3] New TFCSEnergyInterpolationPiecewiseLinear class for piecewise linear interpolation --- .../ISF_FastCaloSimEvent/CMakeLists.txt | 3 +- .../ISF_FastCaloSimEvent/LinkDef.h | 2 + .../TFCSEnergyInterpolationPiecewiseLinear.h | 56 ++++++ ...TFCSEnergyInterpolationPiecewiseLinear.cxx | 169 ++++++++++++++++++ 4 files changed, 229 insertions(+), 1 deletion(-) create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyInterpolationPiecewiseLinear.h create mode 100644 Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyInterpolationPiecewiseLinear.cxx diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt index a0076d6efbe..3a3d33bac99 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt @@ -22,7 +22,7 @@ atlas_depends_on_subdirs( # External dependencies: find_package( CLHEP ) find_package( HepPDT ) -find_package( ROOT COMPONENTS Core Tree MathCore Hist RIO Matrix Physics ) +find_package( ROOT COMPONENTS Core Tree MathCore MathMore Hist RIO Matrix Physics ) find_package( lwtnn REQUIRED ) find_package( LibXml2 ) @@ -55,6 +55,7 @@ atlas_add_root_dictionary( ISF_FastCaloSimEvent _dictSource ISF_FastCaloSimEvent/TFCSInvisibleParametrization.h ISF_FastCaloSimEvent/TFCSInitWithEkin.h ISF_FastCaloSimEvent/TFCSEnergyInterpolationLinear.h + ISF_FastCaloSimEvent/TFCSEnergyInterpolationPiecewiseLinear.h ISF_FastCaloSimEvent/TFCSEnergyInterpolationSpline.h ISF_FastCaloSimEvent/TFCSParametrizationChain.h ISF_FastCaloSimEvent/TFCSParametrizationBinnedChain.h diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/LinkDef.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/LinkDef.h index b9fda2c5725..60d9b266b14 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/LinkDef.h +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/LinkDef.h @@ -27,6 +27,7 @@ #include "ISF_FastCaloSimEvent/TFCSInvisibleParametrization.h" #include "ISF_FastCaloSimEvent/TFCSInitWithEkin.h" #include "ISF_FastCaloSimEvent/TFCSEnergyInterpolationLinear.h" +#include "ISF_FastCaloSimEvent/TFCSEnergyInterpolationPiecewiseLinear.h" #include "ISF_FastCaloSimEvent/TFCSEnergyInterpolationSpline.h" #include "ISF_FastCaloSimEvent/TFCSParametrizationChain.h" #include "ISF_FastCaloSimEvent/TFCSParametrizationBinnedChain.h" @@ -295,6 +296,7 @@ #pragma link C++ class TFCSInvisibleParametrization+; #pragma link C++ class TFCSInitWithEkin+; #pragma link C++ class TFCSEnergyInterpolationLinear+; +#pragma link C++ class TFCSEnergyInterpolationPiecewiseLinear+; #pragma link C++ class TFCSEnergyInterpolationSpline+; #pragma link C++ class TFCSParametrizationChain-; #pragma link C++ class TFCSParametrizationBinnedChain+; diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyInterpolationPiecewiseLinear.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyInterpolationPiecewiseLinear.h new file mode 100644 index 00000000000..24d68c2f2cd --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyInterpolationPiecewiseLinear.h @@ -0,0 +1,56 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef ISF_FASTCALOSIMEVENT_TFCSEnergyInterpolationPiecewiseLinear_h +#define ISF_FASTCALOSIMEVENT_TFCSEnergyInterpolationPiecewiseLinear_h + +//base class include +#include "ISF_FastCaloSimEvent/TFCSParametrization.h" +//interpolator include +#include +//graphic includes for unit_test +#include "TCanvas.h" +#include "TGraph.h" +#include "TAxis.h" + +class TFCSEnergyInterpolationPiecewiseLinear: public TFCSParametrization { + +public: + + TFCSEnergyInterpolationPiecewiseLinear(const char* name = nullptr, const char* title = nullptr); + ~TFCSEnergyInterpolationPiecewiseLinear(); + ///Status bit for energy initialization + enum FCSEnergyInitializationStatusBits { + kOnlyScaleEnergy = BIT(15) ///< Set this bit in the TObject bit field the simulated energy should only be scaled by the spline + }; + + bool OnlyScaleEnergy() const {return TestBit(kOnlyScaleEnergy);}; + void set_OnlyScaleEnergy() {SetBit(kOnlyScaleEnergy);}; + void reset_OnlyScaleEnergy() {ResetBit(kOnlyScaleEnergy);}; + + virtual bool is_match_Ekin_bin(int /*Ekin_bin*/) const override {return true;}; + virtual bool is_match_calosample(int /*calosample*/) const override {return true;}; + + void InitFromArrayInLogEkin(Int_t np, Double_t logEkin[], Double_t response[]); + void InitFromArrayInEkin(Int_t np, Double_t Ekin[], Double_t response[]); + + virtual FCSReturnCode simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const override; + + void Print(Option_t *option="") const override; + + static void unit_test(TFCSSimulationState* simulstate = nullptr,TFCSTruthState* truth = nullptr, const TFCSExtrapolationState* extrapol = nullptr, TGraph* grlinear = nullptr); + +private: + + ROOT::Math::Interpolator * m_linInterpol; + double m_minXValue; + double m_maxXValue; + int m_nPoints; + + ClassDefOverride(TFCSEnergyInterpolationPiecewiseLinear, 1) //TFCSEnergyInterpolationPiecewiseLinear + +}; + +#endif + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyInterpolationPiecewiseLinear.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyInterpolationPiecewiseLinear.cxx new file mode 100644 index 00000000000..8b7993ea2a4 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyInterpolationPiecewiseLinear.cxx @@ -0,0 +1,169 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +*/ + +#include "ISF_FastCaloSimEvent/TFCSEnergyInterpolationPiecewiseLinear.h" +#include "ISF_FastCaloSimEvent/TFCSSimulationState.h" +#include "ISF_FastCaloSimEvent/TFCSTruthState.h" +#include "ISF_FastCaloSimEvent/TFCSExtrapolationState.h" + +#ifdef __FastCaloSimStandAlone__ +namespace Gaudi { + namespace Units { + constexpr double megaelectronvolt = 1.; + constexpr double kiloelectronvolt = 1.e-3 * megaelectronvolt; + constexpr double keV = kiloelectronvolt; + } +} +#else +#include "GaudiKernel/SystemOfUnits.h" +#endif + +//======================================================== +//======= TFCSEnergyInterpolationPiecewiseLinear ========= +//======================================================== + +TFCSEnergyInterpolationPiecewiseLinear::TFCSEnergyInterpolationPiecewiseLinear(const char* name, const char* title): TFCSParametrization(name,title) +{ + ROOT::Math::Interpolator * linInterpol = new ROOT::Math::Interpolator(0, ROOT::Math::Interpolation::kLINEAR); + m_linInterpol = linInterpol; +} + +TFCSEnergyInterpolationPiecewiseLinear::~TFCSEnergyInterpolationPiecewiseLinear() +{ + delete m_linInterpol; +} + +void TFCSEnergyInterpolationPiecewiseLinear::InitFromArrayInLogEkin(Int_t np, Double_t logEkin[], Double_t response[]) +{ + m_linInterpol->SetData(np, logEkin, response); + m_minXValue = *std::min_element(logEkin, logEkin + np); + m_maxXValue = *std::max_element(logEkin, logEkin + np); + m_nPoints = np; +} + +void TFCSEnergyInterpolationPiecewiseLinear::InitFromArrayInEkin(Int_t np, Double_t Ekin[], Double_t response[]) +{ + std::vector logEkin(np); + for(int i=0; iEkin(); + const float Einit = OnlyScaleEnergy() ? simulstate.E() : Ekin; + + //catch very small values of Ekin (use 1 keV here) and fix the interpolation lookup to the 1keV value + const float logEkin = Ekin > Gaudi::Units::keV ? TMath::Log(Ekin) : TMath::Log(Gaudi::Units::keV); + + float Emean; + if(logEkin < m_minXValue){ + Emean = m_linInterpol->Eval(m_minXValue)*Einit; + } + else if(logEkin > m_maxXValue){ + Emean = (m_linInterpol->Eval(m_maxXValue) + m_linInterpol->Deriv(m_maxXValue) * (logEkin - m_maxXValue))*Einit; + } + else{ + Emean = m_linInterpol->Eval(logEkin)*Einit; + } + + if(OnlyScaleEnergy()) ATH_MSG_DEBUG("set E="<=0; + bool longprint = msgLvl(MSG::DEBUG) || (msgLvl(MSG::INFO) && !shortprint); + TString optprint = opt; + optprint.ReplaceAll("short",""); + TFCSParametrization::Print(option); + + if(longprint) ATH_MSG_INFO(optprint <<(OnlyScaleEnergy()?" E()*":" Ekin()*")<<"linInterpol N="<Clone(); + grdraw->SetMarkerColor(46); + grdraw->SetMarkerStyle(8); + + TFCSEnergyInterpolationPiecewiseLinear test("testTFCSEnergyInterpolationPieceWiseLinear", "test TFCSEnergyInterpolationPiecewiseLinear"); + test.set_pdgid(22); + test.set_Ekin_nominal(0.5*(grdraw->GetX()[0]+grdraw->GetX()[grdraw->GetN()-1])); + test.set_Ekin_min(grdraw->GetX()[0]); + test.set_Ekin_max(grdraw->GetX()[grdraw->GetN()-1]); + test.set_eta_nominal(0.225); + test.set_eta_min(0.2); + test.set_eta_max(0.25); + test.InitFromArrayInEkin(grlinear->GetN(),grlinear->GetX(),grlinear->GetY()); + //test.set_OnlyScaleEnergy(); + test.Print(); + + truth->set_pdgid(22); + + TGraph* gr=new TGraph(); + gr->SetNameTitle("testTFCSEnergyInterpolationPiecewiseLogX","test TFCSEnergyInterpolationPiecewiseLinear log x-axis"); + gr->GetXaxis()->SetTitle("Ekin [MeV]"); + gr->GetYaxis()->SetTitle("/Ekin(true)"); + + int ip=0; + for(float Ekin=test.Ekin_min()*0.25;Ekin<=test.Ekin_max()*4;Ekin*=1.05) { + //Init LorentzVector for truth. For photon Ekin=E + truth->SetPxPyPzE(Ekin, 0, 0, Ekin); + simulstate->set_E(Ekin); + if (test.simulate(*simulstate,truth,extrapol) != FCSSuccess) { + return; + } + gr->SetPoint(ip,Ekin,simulstate->E()/Ekin); + ++ip; + } + + //Drawing doesn't make sense inside athena and necessary libraries not linked by default + #if defined(__FastCaloSimStandAlone__) + TCanvas* c = new TCanvas(gr->GetName(),gr->GetTitle()); + gr->Draw("APL"); + grdraw->Draw("Psame"); + c->SetLogx(); + #endif +} \ No newline at end of file -- GitLab From ffa9240b6b1ab215427a15a2ea1d5a635d69c6c4 Mon Sep 17 00:00:00 2001 From: Joshua Falco Beirer Date: Wed, 11 Nov 2020 16:44:22 +0100 Subject: [PATCH 2/3] Made m_linInterpol direct class member --- .../TFCSEnergyInterpolationPiecewiseLinear.h | 3 +-- .../TFCSEnergyInterpolationPiecewiseLinear.cxx | 17 +++++------------ 2 files changed, 6 insertions(+), 14 deletions(-) diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyInterpolationPiecewiseLinear.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyInterpolationPiecewiseLinear.h index 24d68c2f2cd..8bb87367aa2 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyInterpolationPiecewiseLinear.h +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyInterpolationPiecewiseLinear.h @@ -19,7 +19,6 @@ class TFCSEnergyInterpolationPiecewiseLinear: public TFCSParametrization { public: TFCSEnergyInterpolationPiecewiseLinear(const char* name = nullptr, const char* title = nullptr); - ~TFCSEnergyInterpolationPiecewiseLinear(); ///Status bit for energy initialization enum FCSEnergyInitializationStatusBits { kOnlyScaleEnergy = BIT(15) ///< Set this bit in the TObject bit field the simulated energy should only be scaled by the spline @@ -43,7 +42,7 @@ public: private: - ROOT::Math::Interpolator * m_linInterpol; + ROOT::Math::Interpolator m_linInterpol; double m_minXValue; double m_maxXValue; int m_nPoints; diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyInterpolationPiecewiseLinear.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyInterpolationPiecewiseLinear.cxx index 8b7993ea2a4..6a512e9677c 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyInterpolationPiecewiseLinear.cxx +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyInterpolationPiecewiseLinear.cxx @@ -23,20 +23,13 @@ namespace Gaudi { //======= TFCSEnergyInterpolationPiecewiseLinear ========= //======================================================== -TFCSEnergyInterpolationPiecewiseLinear::TFCSEnergyInterpolationPiecewiseLinear(const char* name, const char* title): TFCSParametrization(name,title) +TFCSEnergyInterpolationPiecewiseLinear::TFCSEnergyInterpolationPiecewiseLinear(const char* name, const char* title): TFCSParametrization(name, title), m_linInterpol(0, ROOT::Math::Interpolation::kLINEAR) { - ROOT::Math::Interpolator * linInterpol = new ROOT::Math::Interpolator(0, ROOT::Math::Interpolation::kLINEAR); - m_linInterpol = linInterpol; -} - -TFCSEnergyInterpolationPiecewiseLinear::~TFCSEnergyInterpolationPiecewiseLinear() -{ - delete m_linInterpol; } void TFCSEnergyInterpolationPiecewiseLinear::InitFromArrayInLogEkin(Int_t np, Double_t logEkin[], Double_t response[]) { - m_linInterpol->SetData(np, logEkin, response); + m_linInterpol.SetData(np, logEkin, response); m_minXValue = *std::min_element(logEkin, logEkin + np); m_maxXValue = *std::max_element(logEkin, logEkin + np); m_nPoints = np; @@ -60,13 +53,13 @@ FCSReturnCode TFCSEnergyInterpolationPiecewiseLinear::simulate(TFCSSimulationSta float Emean; if(logEkin < m_minXValue){ - Emean = m_linInterpol->Eval(m_minXValue)*Einit; + Emean = m_linInterpol.Eval(m_minXValue)*Einit; } else if(logEkin > m_maxXValue){ - Emean = (m_linInterpol->Eval(m_maxXValue) + m_linInterpol->Deriv(m_maxXValue) * (logEkin - m_maxXValue))*Einit; + Emean = (m_linInterpol.Eval(m_maxXValue) + m_linInterpol.Deriv(m_maxXValue) * (logEkin - m_maxXValue))*Einit; } else{ - Emean = m_linInterpol->Eval(logEkin)*Einit; + Emean = m_linInterpol.Eval(logEkin)*Einit; } if(OnlyScaleEnergy()) ATH_MSG_DEBUG("set E="<