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 426fdf0c55fd6a724b746829db2fd77bed73b5a4..e65b6f24cccc15ac94b0916d80827a869f87244a 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 c9556bdf6cae84a7e3152fc5899a68b9dfbbd06f..0459e2988b27b14ae011f0a24521339a319e58e6 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 a235190d6fbc1ced7146c607a778e8281d84ca89..2ec25babf3623d0ad27594ab4f884d084f215c9b 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 cfe5a29927dbab5247231500107f698ad4a9d64e..ba3ad8c4ff3d4e27226733eed547d2df2ef6796f 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;
}