diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt index 814e201541956b3b7fb93c94048d5d5fc5cbbb1f..45d8d29cb4d3a77f9620420b0607c40a0b0c71ce 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt +++ b/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 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 03c4fe28ed5656207ddaee6dd82289fd96d41d42..8567ec9c6b3c392056a5e0ca8c0fe29527e42515 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/LinkDef.h +++ b/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+; diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyInterpolationHistogram.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyInterpolationHistogram.h new file mode 100644 index 0000000000000000000000000000000000000000..a730b92e93431ebca8515c8afd0354d1a39a8243 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyInterpolationHistogram.h @@ -0,0 +1,48 @@ +/* + 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 + diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyInterpolationHistogram.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyInterpolationHistogram.cxx new file mode 100644 index 0000000000000000000000000000000000000000..90937730936ae0385aef6916ba9a91534fff2c4a --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyInterpolationHistogram.cxx @@ -0,0 +1,147 @@ +/* + 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 +} +