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 f694f282e8e09798026b160ac4f67e89c4b9aa79..1f9f9d9ce5a723f7e73d7503df5689dfdae643bc 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h
@@ -36,6 +36,9 @@ class TFCSPCAEnergyParametrization:public TFCSEnergyParametrization
void clean();
void Print(Option_t *option = "") const;
+
+ int do_rescale;
+
private:
std::vector<int> m_RelevantLayers;
@@ -45,7 +48,6 @@ class TFCSPCAEnergyParametrization:public TFCSEnergyParametrization
std::vector<TVectorD*> m_SigmaValues;
std::vector<TVectorD*> m_Gauss_means;
std::vector<TVectorD*> m_Gauss_rms;
- std::vector<TVectorD*> m_LowerBounds;
std::vector<std::vector<TFCS1DFunction*> > m_cumulative;
int m_numberpcabins;
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSSimulationState.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSSimulationState.h
index 2aabffd5d981d8fb8534e852afee02a353e18c47..f2ef1a4dc74ad6c7c71f93003bd2b8a6ee8204e4 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSSimulationState.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSSimulationState.h
@@ -36,12 +36,15 @@ class TFCSSimulationState:public TObject
void deposit(const CaloDetDescrElement* cellele, float E);
void Print(Option_t *option="") const;
+ void set_SF(double mysf) {SF=mysf;}
+ double get_SF() {return SF;}
void clear();
private:
int m_Ebin;
double m_Etot;
// TO BE CLEANED UP! SHOULD ONLY STORE EITHER E OR EFRAC!!!
+ double SF;
double m_E[CaloCell_ID_FCS::MaxSample];
double m_Efrac[CaloCell_ID_FCS::MaxSample];
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
index 3007f1c27cd9f95bbd846e60c78515e883f74689..bf3695b474dcfc300ea31d6a24d7215518233c31 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
@@ -23,6 +23,7 @@
TFCSPCAEnergyParametrization::TFCSPCAEnergyParametrization(const char* name, const char* title):TFCSEnergyParametrization(name,title)
{
m_numberpcabins=1;
+ do_rescale=1;
}
bool TFCSPCAEnergyParametrization::is_match_Ekin_bin(int Ekin_bin) const
@@ -67,97 +68,71 @@ void TFCSPCAEnergyParametrization::simulate(TFCSSimulationState& simulstate,cons
TVectorD* SigmaValues =m_SigmaValues[pcabin-1];
TVectorD* Gauss_means =m_Gauss_means[pcabin-1];
TVectorD* Gauss_rms =m_Gauss_rms[pcabin-1];
- TVectorD* LowerBounds =m_LowerBounds[pcabin-1];
std::vector<TFCS1DFunction*> cumulative=m_cumulative[pcabin-1];
TRandom3* Random3=new TRandom3(); Random3->SetSeed(0);
- std::vector<std::string> layer;
std::vector<int> layerNr;
for(unsigned int i=0;i<m_RelevantLayers.size();i++)
- layerNr.push_back(m_RelevantLayers[i]);
- for(unsigned int i=0;i<layerNr.size();i++)
- {
- std::string thislayer=Form("layer%i",layerNr[i]);
- layer.push_back(thislayer);
- }
- layer.push_back("totalE");
+ layerNr.push_back(m_RelevantLayers[i]);
double* vals_gauss_means=(double*)Gauss_means->GetMatrixArray();
double* vals_gauss_rms =Gauss_rms->GetMatrixArray();
- double* vals_lowerBounds=LowerBounds->GetMatrixArray();
- double *output_data = new double[layer.size()] ;
- double *input_data = new double[layer.size()] ;
+ double *output_data = new double[layerNr.size()+1];
+ double *input_data = new double[layerNr.size()+1];
- for(unsigned int l=0;l<layer.size();l++)
- {
- double mean=vals_gauss_means[l];
- double rms =vals_gauss_rms[l];
- double gauszz=Random3->Gaus(mean,rms);
- input_data[l]=gauszz;
- }
+ for(unsigned int l=0;l<=layerNr.size();l++)
+ {
+ double mean=vals_gauss_means[l];
+ double rms =vals_gauss_rms[l];
+ double gauszz=Random3->Gaus(mean,rms);
+ input_data[l]=gauszz;
+ }
- P2X(SigmaValues, MeanValues, EV, layer.size(), input_data, output_data, layer.size());
+ P2X(SigmaValues, MeanValues, EV, layerNr.size()+1, input_data, output_data, layerNr.size()+1);
- double *simdata_uniform = new double[layer.size()];
- double *simdata = new double[layer.size()];
- double *simdata_scaled = new double[layer.size()];
+ double *simdata = new double[layerNr.size()+1];
double sum_fraction=0.0;
- for(unsigned int l=0;l<layer.size();l++)
- {
- simdata_uniform[l]=(TMath::Erf(output_data[l]/1.414213562)+1)/2.f;
-
- if(simdata_uniform[l]<vals_lowerBounds[l]) simdata_uniform[l]=vals_lowerBounds[l];
-
- simdata[l]=cumulative[l]->rnd_to_fct(simdata_uniform[l]);
-
- if(simdata[l]<0) simdata[l]=0;
- if(layer[l]!="totalE" && simdata[l]>1) simdata[l]=1;
- if(layer[l]!="totalE") sum_fraction+=simdata[l];
-
- }
-
- double scale=1.0/sum_fraction;
-
- sum_fraction=0.0;
-
- for(unsigned int l=0;l<layer.size();l++)
- {
- if(layer[l]!="totalE")
- {
- simdata_scaled[l]=simdata[l]*scale;
- if(l<layerNr.size())
- {
- sum_fraction+=simdata_scaled[l];
- }
- }
- }
-
- double total_energy=simdata[layer.size()-1]*simulstate.E()/Ekin_nominal();
- //double total_energy=simdata[layer.size()-1];
+ 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);
+
+ 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);
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_scaled[l];
- }
- simulstate.set_Efrac(s,energyfrac);
- simulstate.set_E(s,energyfrac*total_energy);
- }
+ {
+ 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);
+ }
delete Random3;
delete [] output_data;
delete [] input_data;
delete [] simdata;
- delete [] simdata_uniform;
- delete [] simdata_scaled;
-
}
@@ -228,14 +203,12 @@ bool TFCSPCAEnergyParametrization::loadInputs(TFile* file, std::string folder)
TVectorD* SigmaValues =(TVectorD*)gDirectory->Get("SigmaValues");
TVectorD* Gauss_means =(TVectorD*)gDirectory->Get("Gauss_means");
TVectorD* Gauss_rms =(TVectorD*)gDirectory->Get("Gauss_rms");
- TVectorD* LowerBounds =(TVectorD*)gDirectory->Get("LowerBounds");
if(symCov == NULL) {ATH_MSG_WARNING("TFCSPCAEnergyParametrization::symCov in pcabin "<<bin<<" is null!"); load_ok=false;}
if(MeanValues == NULL) {ATH_MSG_WARNING("TFCSPCAEnergyParametrization::MeanValues in pcabin "<<bin<<" is null!"); load_ok=false;}
if(SigmaValues == NULL) {ATH_MSG_WARNING("TFCSPCAEnergyParametrization::SigmaValues in pcabin "<<bin<<" is null!"); load_ok=false;}
if(Gauss_means == NULL) {ATH_MSG_WARNING("TFCSPCAEnergyParametrization::Gauss_means in pcabin "<<bin<<" is null!"); load_ok=false;}
if(Gauss_rms == NULL) {ATH_MSG_WARNING("TFCSPCAEnergyParametrization::Gause_rms in pcabin "<<bin<<" is null!"); load_ok=false;}
- if(LowerBounds == NULL) {ATH_MSG_WARNING("TFCSPCAEnergyParametrization::LowerBounds in pcabin "<<bin<<" is null!"); load_ok=false;}
if(!load_ok) return false;
@@ -246,7 +219,6 @@ bool TFCSPCAEnergyParametrization::loadInputs(TFile* file, std::string folder)
m_SigmaValues.push_back(SigmaValues);
m_Gauss_means.push_back(Gauss_means);
m_Gauss_rms.push_back(Gauss_rms);
- m_LowerBounds.push_back(LowerBounds);
std::vector<std::string> layer;
std::vector<int> layerNr;