Merge branch 'cherry-pick-f96b9017903-master' into 'master'

Sweeping !39061 from 21.0 to master.
21.0 add tfcs energy interpolation histogram

See merge request atlas/athena!39337

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt

@@ -41,6 +41,7 @@ atlas_add_root_dictionary( ISF_FastCaloSimEvent _dictSource
                            ISF_FastCaloSimEvent/TFCSParametrization.h 
                            ISF_FastCaloSimEvent/TFCSInvisibleParametrization.h
                            ISF_FastCaloSimEvent/TFCSInitWithEkin.h
+                           ISF_FastCaloSimEvent/TFCSEnergyInterpolationHistogram.h
                            ISF_FastCaloSimEvent/TFCSEnergyInterpolationLinear.h
                            ISF_FastCaloSimEvent/TFCSEnergyInterpolationPiecewiseLinear.h
                            ISF_FastCaloSimEvent/TFCSEnergyInterpolationSpline.h

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/LinkDef.h

@@ -26,6 +26,7 @@
 #include "ISF_FastCaloSimEvent/TFCSParametrization.h"
 #include "ISF_FastCaloSimEvent/TFCSInvisibleParametrization.h"
 #include "ISF_FastCaloSimEvent/TFCSInitWithEkin.h"
+#include "ISF_FastCaloSimEvent/TFCSEnergyInterpolationHistogram.h"
 #include "ISF_FastCaloSimEvent/TFCSEnergyInterpolationLinear.h"
 #include "ISF_FastCaloSimEvent/TFCSEnergyInterpolationPiecewiseLinear.h"
 #include "ISF_FastCaloSimEvent/TFCSEnergyInterpolationSpline.h"
@@ -295,6 +296,7 @@
 #pragma link C++ class TFCSParametrization+;
 #pragma link C++ class TFCSInvisibleParametrization+;
 #pragma link C++ class TFCSInitWithEkin+;
+#pragma link C++ class TFCSEnergyInterpolationHistogram+;
 #pragma link C++ class TFCSEnergyInterpolationLinear+;
 #pragma link C++ class TFCSEnergyInterpolationPiecewiseLinear-;
 #pragma link C++ class TFCSEnergyInterpolationSpline+;

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyInterpolationHistogram.h

+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef ISF_FASTCALOSIMEVENT_TFCSEnergyInterpolationHistogram_h
+#define ISF_FASTCALOSIMEVENT_TFCSEnergyInterpolationHistogram_h
+
+#include "ISF_FastCaloSimEvent/TFCSParametrization.h"
+#include "TH1F.h"
+
+class TFCSEnergyInterpolationHistogram:public TFCSParametrization {
+public:
+  TFCSEnergyInterpolationHistogram(const char* name=nullptr, const char* title=nullptr);
+
+  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;};
+  
+  ///Initialize interpolation from histogram
+  ///x values should be Ekin, y values should <E(reco)/Ekin(true)>
+  void InitFromHist(const TH1F& hist) {m_hist=hist;};
+
+  const TH1F& hist() const {return m_hist;};
+
+  ///Initialize simulstate with the mean reconstructed energy in the calorimater expeted from the true kinetic energy
+  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,TH1F* hist=nullptr);
+private:
+  TH1F m_hist;
+
+  ClassDefOverride(TFCSEnergyInterpolationHistogram,1)  //TFCSEnergyInterpolationHistogram
+};
+
+#if defined(__ROOTCLING__) && defined(__FastCaloSimStandAlone__)
+#pragma link C++ class TFCSEnergyInterpolationHistogram+;
+#endif
+
+#endif
+

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

+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "ISF_FastCaloSimEvent/TFCSEnergyInterpolationHistogram.h"
+#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
+#include "ISF_FastCaloSimEvent/TFCSTruthState.h"
+#include "ISF_FastCaloSimEvent/TFCSExtrapolationState.h"
+#include "TFile.h"
+#include "TCanvas.h"
+#include "TGraph.h"
+#include "TAxis.h"
+#include <iostream>
+#include <vector>
+
+#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
+
+//=============================================
+//======= TFCSEnergyInterpolation =========
+//=============================================
+
+TFCSEnergyInterpolationHistogram::TFCSEnergyInterpolationHistogram(const char* name, const char* title):TFCSParametrization(name,title)
+{
+}
+
+FCSReturnCode TFCSEnergyInterpolationHistogram::simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState*) const
+{
+  float Emean;
+  float Einit;
+  const float Ekin=truth->Ekin();
+  if(OnlyScaleEnergy()) Einit=simulstate.E();
+  else Einit=Ekin;
+  if(Einit<m_hist.GetXaxis()->GetBinLowEdge(1)) {
+    Emean=m_hist.GetBinContent(1)*Einit;
+  } else {
+    if(Einit>m_hist.GetXaxis()->GetBinUpEdge(m_hist.GetNbinsX())) {
+      Emean=m_hist.GetBinContent(m_hist.GetNbinsX())*Einit;
+    } else {
+      Emean=m_hist.GetBinContent(m_hist.GetXaxis()->FindBin(Einit))*Einit;
+    }  
+  }  
+
+  if(OnlyScaleEnergy()) {
+    ATH_MSG_DEBUG("set E="<<Emean<<" for true Ekin="<<truth->Ekin()<<" and E="<<Einit);
+  }
+  else{
+    ATH_MSG_DEBUG("set E="<<Emean<<" for true Ekin="<<truth->Ekin());
+  }
+  simulstate.set_E(Emean);
+
+  return FCSSuccess;
+}
+
+void TFCSEnergyInterpolationHistogram::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()*")<<"histNbins="<<m_hist.GetNbinsX()
+                           <<" "<<m_hist.GetXaxis()->GetBinLowEdge(1)<<"<=Ekin<="<<m_hist.GetXaxis()->GetBinUpEdge(m_hist.GetNbinsX()));
+}
+
+void TFCSEnergyInterpolationHistogram::unit_test(TFCSSimulationState* simulstate,TFCSTruthState* truth, const TFCSExtrapolationState* extrapol,TH1F* hist)
+{
+  if(!simulstate) simulstate=new TFCSSimulationState();
+  if(!truth) truth=new TFCSTruthState();
+  if(!extrapol) extrapol=new TFCSExtrapolationState();
+  
+  if(!hist) {
+    hist = new TH1F("h1", "h1 title", 16, 0., 512.);
+    hist->SetBinContent( 1,0.687165);
+    hist->SetBinContent( 2,0.756837);
+    hist->SetBinContent( 3,0.836673);
+    hist->SetBinContent( 4,0.896336);
+    hist->SetBinContent( 5,0.889785);
+    hist->SetBinContent( 6,0.901266);
+    hist->SetBinContent( 7,0.888937);
+    hist->SetBinContent( 8,0.919943);
+    hist->SetBinContent( 9,0.941806);
+    hist->SetBinContent(10,0.934668);
+    hist->SetBinContent(11,0.939502);
+    hist->SetBinContent(12,0.940796);
+    hist->SetBinContent(13,0.945894);
+    hist->SetBinContent(14,0.955445);
+    hist->SetBinContent(15,0.955593);
+    hist->SetBinContent(16,0.943673);
+  }  
+  
+  TH1F* hdraw=(TH1F*)hist->Clone();
+  hdraw->SetMarkerColor(46);
+  hdraw->SetMarkerStyle(8);
+  
+  TFCSEnergyInterpolationHistogram test("testTFCSEnergyInterpolationHistogram","test TFCSEnergyInterpolationHistogram");
+  test.set_pdgid(22);
+  test.set_Ekin_nominal(0.5*(hdraw->GetXaxis()->GetBinLowEdge(1)+hdraw->GetXaxis()->GetBinUpEdge(hdraw->GetNbinsX())));
+  test.set_Ekin_min(hdraw->GetXaxis()->GetBinLowEdge(1));
+  test.set_Ekin_max(hdraw->GetXaxis()->GetBinUpEdge(hdraw->GetNbinsX()));
+  test.set_eta_nominal(0.225);
+  test.set_eta_min(0.2);
+  test.set_eta_max(0.25);
+  test.InitFromHist(*hist);
+  //test.set_OnlyScaleEnergy();
+  test.Print();
+  test.hist().Dump();
+  
+  truth->set_pdgid(22);
+  
+  TGraph* gr=new TGraph();
+  gr->SetNameTitle("testTFCSEnergyInterpolationHistogramLogX","test TFCSEnergyInterpolationHistogram");
+  gr->GetXaxis()->SetTitle("Ekin [MeV]");
+  gr->GetYaxis()->SetTitle("<E(reco)>/Ekin(true)");
+  
+  int ip=0;
+  for(float Ekin=std::max(test.Ekin_min()*0.25,0.1);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(hdraw->GetName(),hdraw->GetTitle());
+  hdraw->Draw("HIST");
+  gr->Draw("same,APL");
+  c->SetLogx();
+  #endif
+}
+