Commit f96b9017 authored by Rachid Mazini's avatar Rachid Mazini
Browse files

Merge branch '21.0-addTFCSEnergyInterpolationHistogram' into '21.0'

21.0 add tfcs energy interpolation histogram

See merge request !39061
parents 02b107d3 958d7091
......@@ -54,6 +54,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
......
......@@ -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+;
......
/*
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
/*
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
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment