Merge branch '21.0-PiecewiseLinearEnergyInterpolation' into '21.0'

New TFCSEnergyInterpolationPiecewiseLinear class for piecewise linear interpolation

See merge request atlas/athena!38095

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

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+;

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyInterpolationPiecewiseLinear.h

+/*
+  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 <Math/Interpolator.h>
+//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);
+  ///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, const Double_t logEkin[], const Double_t response[]);
+  void InitFromArrayInEkin(Int_t np, const Double_t Ekin[], const 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
+

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyInterpolationPiecewiseLinear.cxx

+/*
+  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), m_linInterpol(0, ROOT::Math::Interpolation::kLINEAR)
+{
+}
+
+void TFCSEnergyInterpolationPiecewiseLinear::InitFromArrayInLogEkin(Int_t np, const Double_t logEkin[], const 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, const Double_t Ekin[], const Double_t response[])
+{
+  std::vector<Double_t> logEkin(np);
+  for(int i=0; i<np; i++) logEkin[i] = TMath::Log(Ekin[i]);
+  InitFromArrayInLogEkin(np, logEkin.data(), response);
+}
+
+FCSReturnCode TFCSEnergyInterpolationPiecewiseLinear::simulate(TFCSSimulationState& simulstate, const TFCSTruthState* truth, const TFCSExtrapolationState*) const
+{
+
+  const float Ekin  = truth->Ekin();
+  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="<<Emean<<" for true Ekin="<<Ekin<<" and E="<<Einit);
+  else                  ATH_MSG_DEBUG("set E="<<Emean<<" for true Ekin="<<Ekin);
+
+  //set mean energy of simulstate
+  simulstate.set_E(Emean);
+
+  return FCSSuccess;
+}
+
+void TFCSEnergyInterpolationPiecewiseLinear::Print(Option_t *option) const
+{
+  TString opt(option);
+  bool shortprint = opt.Index("short")>=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="<<m_nPoints
+                           <<" "<<m_minXValue<<"<=log(Ekin)<="<<m_maxXValue
+                           <<" "<<TMath::Exp(m_minXValue)<<"<=Ekin<="<<TMath::Exp(m_maxXValue));
+}
+
+
+void TFCSEnergyInterpolationPiecewiseLinear::unit_test(TFCSSimulationState* simulstate, TFCSTruthState* truth, const TFCSExtrapolationState* extrapol, TGraph* grlinear)
+{
+  if(!simulstate) simulstate  = new TFCSSimulationState();
+  if(!truth)      truth       = new TFCSTruthState();
+  if(!extrapol)   extrapol    = new TFCSExtrapolationState();
+
+  if(!grlinear) {
+    const int Graph0_n=9;
+    Double_t Graph0_fx1001[Graph0_n] = {1.024,
+                                        2.048,
+                                        4.094,
+                                        8.192,
+                                        16.384,
+                                        32.768,
+                                        65.536,
+                                        131.072,
+                                        262.144};
+    
+    for(int i=0; i<Graph0_n; ++i) Graph0_fx1001[i]*=1000;
+    Double_t Graph0_fy1001[Graph0_n] = {0.6535402,
+                                        0.6571529,
+                                        0.6843001,
+                                        0.7172835,
+                                        0.7708416,
+                                        0.798819,
+                                        0.8187628,
+                                        0.8332745,
+                                        0.8443931};
+    grlinear = new TGraph(Graph0_n,Graph0_fx1001,Graph0_fy1001);
+  }  
+
+  TGraph* grdraw=(TGraph*)grlinear->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("<E(reco)>/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