Merge branch '21.0-FCS-PCAbinzero-support' into '21.0'

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

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;
 

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;

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)) {

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;
 }