Merge branch 'cherry-pick-7544f777c4-master' into 'master'

Sweeping !15227 from 21.0 to master.
FCS support for PCA bin 0

See merge request atlas/athena!15475

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