From f5864db77b23b9af15fbe9de74b99af96df5b10e Mon Sep 17 00:00:00 2001
From: James Beacham <j.beacham@cern.ch>
Date: Tue, 30 Oct 2018 19:10:10 +0000
Subject: [PATCH] Merge branch '21.0-FCS-PCAbinzero-support' into '21.0'
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
FCS support for PCA bin 0
See merge request atlas/athena!15227
(cherry picked from commit 7544f777c47168daea476f67cfb6e528d5b0a934)
0c31d2e2 PCA bin may now be 0, in case of 0 return zero energies
a73ed150 Adding support to read pca bin probabilities from the second pca file
ca88ea00 PCA bin probs are now loaded in TFCSEnergyBinParametrization
e54a59a3 deleting prob vector from PCAEnergyParametrization and adding size check in…
b11111e1 Option to set PCA bin probabilities in case no pcabinprob object is available
---
.../python/iconfTool/gui/__init__.py | 0
.../TFCSEnergyBinParametrization.h | 2 +
.../TFCSPCAEnergyParametrization.h | 1 +
.../src/TFCSEnergyBinParametrization.cxx | 53 +++++++-
.../src/TFCSPCAEnergyParametrization.cxx | 127 ++++++++++--------
5 files changed, 125 insertions(+), 58 deletions(-)
mode change 100644 => 100755 Control/AthenaConfiguration/python/iconfTool/gui/__init__.py
diff --git a/Control/AthenaConfiguration/python/iconfTool/gui/__init__.py b/Control/AthenaConfiguration/python/iconfTool/gui/__init__.py
old mode 100644
new mode 100755
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyBinParametrization.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyBinParametrization.h
index 426fdf0c55f..e65b6f24ccc 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyBinParametrization.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyBinParametrization.h
@@ -6,6 +6,7 @@
#define ISF_FASTCALOSIMEVENT_TFCSEnergyBinParametrization_h
#include "ISF_FastCaloSimEvent/TFCSEnergyParametrization.h"
+#include "TFile.h"
#include<map>
#include<vector>
@@ -30,6 +31,7 @@ class TFCSEnergyBinParametrization:public TFCSEnergyParametrization
/// the function will normalize probabilities automatically, if the sum of values is not 1
/// current convention is to start Ekin_bin counting at 1, to be updated to start counting with 0
virtual void set_pdgid_Ekin_bin_probability(int id,std::vector< float > prob);
+ virtual void load_pdgid_Ekin_bin_probability_from_file(int id, TFile* file, std::string prob_object_name);
virtual FCSReturnCode simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) override;
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h
index c9556bdf6ca..0459e2988b2 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h
@@ -33,6 +33,7 @@ class TFCSPCAEnergyParametrization:public TFCSEnergyParametrization
void P2X(TVectorD*, TVectorD* , TMatrixD* , int, double* , double* , int);
bool loadInputs(TFile* file);
bool loadInputs(TFile* file,std::string);
+
void clean();
void Print(Option_t *option = "") const override;
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyBinParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyBinParametrization.cxx
index a235190d6fb..2ec25babf36 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyBinParametrization.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyBinParametrization.cxx
@@ -7,7 +7,8 @@
#include "ISF_FastCaloSimEvent/TFCSEnergyBinParametrization.h"
#include "ISF_FastCaloSimEvent/TFCSTruthState.h"
#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
-
+#include "TFile.h"
+#include "TVectorF.h"
#include "TMath.h"
//=============================================
@@ -78,7 +79,55 @@ void TFCSEnergyBinParametrization::set_pdgid_Ekin_bin_probability(int id,std::ve
m_pdgid_Ebin_probability[id][iEbin]=p;
}
}
+
+
+void TFCSEnergyBinParametrization::load_pdgid_Ekin_bin_probability_from_file(int id, TFile* file, std::string prob_object_name)
+{
+ add_pdgid(id);
+
+ float* prob;
+ file->cd();
+ TVectorF* pcabinprobvector=(TVectorF*)gDirectory->Get(prob_object_name.c_str());
+
+ if(!pcabinprobvector)
+ {
+ ATH_MSG_INFO("TFCSEnergyBinParametrization::"<<prob_object_name<<" is null. Using equal PCA probabilities.");
+
+ prob=new float[m_pdgid_Ebin_probability[id].size()];
+
+ prob[0]=0.0;
+
+ for(unsigned int i=1;i<m_pdgid_Ebin_probability[id].size();i++)
+ {
+ prob[i]=1.0;
+ }
+
+ }
+
+ if(pcabinprobvector)
+ {
+ if((unsigned int)pcabinprobvector->GetNoElements()!=m_pdgid_Ebin_probability[id].size())
+ {
+ ATH_MSG_ERROR("TFCSEnergyBinParametrization::set_pdgid_Ekin_bin_probability(): size of prob array does not match! in.size()="<<pcabinprobvector->GetNoElements()<<" instance="<<m_pdgid_Ebin_probability[id].size());
+ return;
+ }
+
+ prob=pcabinprobvector->GetMatrixArray();
+
+ }
+
+ float ptot=0;
+ for(int iEbin=0;iEbin<=n_bins();++iEbin) ptot+=prob[iEbin];
+ float p=0;
+ for(int iEbin=0;iEbin<=n_bins();++iEbin)
+ {
+ p+=prob[iEbin]/ptot;
+ m_pdgid_Ebin_probability[id][iEbin]=p;
+ }
+
+}
+
void TFCSEnergyBinParametrization::Print(Option_t *option) const
{
TString opt(option);
@@ -119,7 +168,7 @@ FCSReturnCode TFCSEnergyBinParametrization::simulate(TFCSSimulationState& simuls
float searchRand = CLHEP::RandFlat::shoot(simulstate.randomEngine());
int chosenBin=TMath::BinarySearch(n_bins()+1, m_pdgid_Ebin_probability[pdgid].data(), searchRand)+1;
- if(chosenBin<1 || chosenBin>n_bins()) {
+ if(chosenBin<0 || chosenBin>n_bins()) {
ATH_MSG_ERROR("TFCSEnergyBinParametrization::simulate(): cannot simulate bin="<<chosenBin);
ATH_MSG_ERROR(" This error could probably be retried.");
if(msgLvl(MSG::ERROR)) {
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
index cfe5a29927d..ba3ad8c4ff3 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
@@ -61,79 +61,94 @@ void TFCSPCAEnergyParametrization::Print(Option_t *option) const
FCSReturnCode TFCSPCAEnergyParametrization::simulate(TFCSSimulationState& simulstate,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* /*extrapol*/)
{
+
if (!simulstate.randomEngine()) {
return FCSFatal;
}
int pcabin=simulstate.Ebin();
-
- TMatrixD* EV =m_EV[pcabin-1];
- TVectorD* MeanValues =m_MeanValues[pcabin-1];
- TVectorD* SigmaValues =m_SigmaValues[pcabin-1];
- TVectorD* Gauss_means =m_Gauss_means[pcabin-1];
- TVectorD* Gauss_rms =m_Gauss_rms[pcabin-1];
- std::vector<TFCS1DFunction*> cumulative=m_cumulative[pcabin-1];
-
- std::vector<int> layerNr;
- for(unsigned int i=0;i<m_RelevantLayers.size();i++)
- layerNr.push_back(m_RelevantLayers[i]);
- double* vals_gauss_means=(double*)Gauss_means->GetMatrixArray();
- double* vals_gauss_rms =Gauss_rms->GetMatrixArray();
-
- double *output_data = new double[layerNr.size()+1];
- double *input_data = new double[layerNr.size()+1];
-
- for(unsigned int l=0;l<=layerNr.size();l++)
+ if(pcabin==0)
{
- double mean=vals_gauss_means[l];
- double rms =vals_gauss_rms[l];
- double gauszz = CLHEP::RandGauss::shoot(simulstate.randomEngine(), mean, rms);
- input_data[l]=gauszz;
+ simulstate.set_E(0);
+ for(int s=0;s<CaloCell_ID_FCS::MaxSample;s++)
+ {
+ simulstate.set_E(s,0.0);
+ simulstate.set_Efrac(s,0.0);
+ }
}
+ else
+ {
+
+ TMatrixD* EV =m_EV[pcabin-1];
+ TVectorD* MeanValues =m_MeanValues[pcabin-1];
+ TVectorD* SigmaValues =m_SigmaValues[pcabin-1];
+ TVectorD* Gauss_means =m_Gauss_means[pcabin-1];
+ TVectorD* Gauss_rms =m_Gauss_rms[pcabin-1];
+ std::vector<TFCS1DFunction*> cumulative=m_cumulative[pcabin-1];
+
+ std::vector<int> layerNr;
+ for(unsigned int i=0;i<m_RelevantLayers.size();i++)
+ layerNr.push_back(m_RelevantLayers[i]);
+
+ double* vals_gauss_means=(double*)Gauss_means->GetMatrixArray();
+ double* vals_gauss_rms =Gauss_rms->GetMatrixArray();
- P2X(SigmaValues, MeanValues, EV, layerNr.size()+1, input_data, output_data, layerNr.size()+1);
+ double *output_data = new double[layerNr.size()+1];
+ double *input_data = new double[layerNr.size()+1];
- double *simdata = new double[layerNr.size()+1];
- double sum_fraction=0.0;
- for(unsigned int l=0;l<=layerNr.size();l++)
- {
- double simdata_uniform=(TMath::Erf(output_data[l]/1.414213562)+1)/2.f;
+ for(unsigned int l=0;l<=layerNr.size();l++)
+ {
+ double mean=vals_gauss_means[l];
+ double rms =vals_gauss_rms[l];
+ double gauszz = CLHEP::RandGauss::shoot(simulstate.randomEngine(), mean, rms);
+ input_data[l]=gauszz;
+ }
+
+ P2X(SigmaValues, MeanValues, EV, layerNr.size()+1, input_data, output_data, layerNr.size()+1);
+
+ double *simdata = new double[layerNr.size()+1];
+ double sum_fraction=0.0;
+ for(unsigned int l=0;l<=layerNr.size();l++)
+ {
+ double simdata_uniform=(TMath::Erf(output_data[l]/1.414213562)+1)/2.f;
+
+ simdata[l]=cumulative[l]->rnd_to_fct(simdata_uniform);
- simdata[l]=cumulative[l]->rnd_to_fct(simdata_uniform);
+ if(l!=layerNr.size()) //sum up the fractions, but not the totalE
+ sum_fraction+=simdata[l];
+ }
- if(l!=layerNr.size()) //sum up the fractions, but not the totalE
- sum_fraction+=simdata[l];
- }
-
- double scalefactor=1.0/sum_fraction;
- if(!do_rescale) scalefactor=1.0;
-
- for(unsigned int l=0;l<layerNr.size();l++)
- {
- simdata[l]*=scalefactor;
- }
-
- double total_energy=simdata[layerNr.size()]*simulstate.E()/Ekin_nominal();
- simulstate.set_E(total_energy);
- ATH_MSG_DEBUG("set E to total_energy="<<total_energy);
+ double scalefactor=1.0/sum_fraction;
+ if(!do_rescale) scalefactor=1.0;
- for(int s=0;s<CaloCell_ID_FCS::MaxSample;s++)
- {
- double energyfrac=0.0;
for(unsigned int l=0;l<layerNr.size();l++)
{
- if(layerNr[l]==s)
- energyfrac=simdata[l];
+ simdata[l]*=scalefactor;
+ }
+
+ double total_energy=simdata[layerNr.size()]*simulstate.E()/Ekin_nominal();
+ simulstate.set_E(total_energy);
+ ATH_MSG_DEBUG("set E to total_energy="<<total_energy);
+
+ for(int s=0;s<CaloCell_ID_FCS::MaxSample;s++)
+ {
+ double energyfrac=0.0;
+ for(unsigned int l=0;l<layerNr.size();l++)
+ {
+ if(layerNr[l]==s)
+ energyfrac=simdata[l];
+ }
+ simulstate.set_Efrac(s,energyfrac);
+ simulstate.set_E(s,energyfrac*total_energy);
+ simulstate.set_SF(scalefactor);
}
- simulstate.set_Efrac(s,energyfrac);
- simulstate.set_E(s,energyfrac*total_energy);
- simulstate.set_SF(scalefactor);
- }
- delete [] output_data;
- delete [] input_data;
- delete [] simdata;
+ delete [] output_data;
+ delete [] input_data;
+ delete [] simdata;
+
+ }
return FCSSuccess;
}
--
GitLab