From 200ebf4cdd980ee591e524ed78335f3aa1693ad4 Mon Sep 17 00:00:00 2001
From: Jana Schaarschmidt <jana.schaarschmidt@cern.ch>
Date: Thu, 19 Apr 2018 20:10:03 +0200
Subject: [PATCH] Updates of FCS energy para and sim
Former-commit-id: c79a7ed1b54b8c6125823d5437431afb2c67a0fb
---
.../TFCS1DFunctionHistogram.h | 2 +-
.../TFCS1DFunctionRegressionTF.h | 2 +
.../TFCSPCAEnergyParametrization.h | 6 +-
.../src/TFCS1DFunctionHistogram.cxx | 60 ++++++++------
.../src/TFCSPCAEnergyParametrization.cxx | 82 +++++++++++--------
5 files changed, 87 insertions(+), 65 deletions(-)
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunctionHistogram.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunctionHistogram.h
index dc052a790fb8..16310ed76c18 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunctionHistogram.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunctionHistogram.h
@@ -23,6 +23,7 @@ class TFCS1DFunctionHistogram:public TFCS1DFunction
TH1* vector_to_histo();
double get_inverse(double rnd);
double linear(double x1,double x2,double y1,double y2,double x);
+ double non_linear(double x1,double x2,double y1,double y2,double x);
double get_maxdev(TH1*, TH1D*);
void smart_rebin_loop(TH1* hist, double);
@@ -40,7 +41,6 @@ class TFCS1DFunctionHistogram:public TFCS1DFunction
vector<float> m_HistoBorders;
vector<float> m_HistoContents;
-
ClassDef(TFCS1DFunctionHistogram,1) //TFCS1DFunctionHistogram
};
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunctionRegressionTF.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunctionRegressionTF.h
index e3d9d2aa0d1f..34e387d83a00 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunctionRegressionTF.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunctionRegressionTF.h
@@ -31,9 +31,11 @@ class TFCS1DFunctionRegressionTF:public TFCS1DFunctionRegression
ClassDef(TFCS1DFunctionRegressionTF,1)
};
+
#if defined(__ROOTCLING__) && defined(__FastCaloSimStandAlone__)
#pragma link C++ class TFCS1DFunctionRegressionTF+;
#endif
#endif
+
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 9b0d01007eae..34d6c1acc620 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h
@@ -31,15 +31,15 @@ class TFCSPCAEnergyParametrization:public TFCSEnergyParametrization
virtual bool is_match_all_calosample() const {return false;};
void P2X(TVectorD*, TVectorD* , TMatrixD* , int, double* , double* , int);
- void loadInputs(TFile* file);
- void loadInputs(TFile* file,std::string);
+ bool loadInputs(TFile* file);
+ bool loadInputs(TFile* file,std::string);
void Print(Option_t *option = "") const;
private:
std::vector<int> m_RelevantLayers;
- std::vector<TMatrixDSym*> m_symCov;
+ std::vector<TMatrixD*> m_EV;
std::vector<TVectorD*> m_MeanValues;
std::vector<TVectorD*> m_SigmaValues;
std::vector<TVectorD*> m_Gauss_means;
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCS1DFunctionHistogram.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCS1DFunctionHistogram.cxx
index 9294c8ccebb2..70846c31a0d9 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCS1DFunctionHistogram.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCS1DFunctionHistogram.cxx
@@ -113,7 +113,7 @@ void TFCS1DFunctionHistogram::smart_rebin_loop(TH1* hist, double cut_maxdev)
maxdev=get_maxdev(hist,h_out);
maxdev*=100.0;
-
+
if(i%100==0) cout<<"Iteration nr. "<<i<<" -----> change "<<change<<" bins "<<h_out->GetNbinsX()<<" -> maxdev="<<maxdev<<endl;
if(maxdev<cut_maxdev && h_out->GetNbinsX()>5 && i<1000)
@@ -136,14 +136,14 @@ void TFCS1DFunctionHistogram::smart_rebin_loop(TH1* hist, double cut_maxdev)
cout<<"Info: Rebinned histogram has "<<h_output->GetNbinsX()<<" bins."<<endl;
//store:
- for(int b=1;b<=h_output->GetNbinsX();b++)
- m_HistoContents.push_back(h_output->GetBinContent(b));
for(int b=1;b<=h_output->GetNbinsX();b++)
m_HistoBorders.push_back((float)h_output->GetBinLowEdge(b));
-
m_HistoBorders.push_back((float)h_output->GetXaxis()->GetXmax());
+ for(int b=1;b<h_output->GetNbinsX();b++)
+ m_HistoContents.push_back(h_output->GetBinContent(b));
+ m_HistoContents.push_back(1);
}
@@ -244,6 +244,7 @@ double TFCS1DFunctionHistogram::rnd_to_fct(double rnd)
{
double value2=get_inverse(rnd);
+
return value2;
}
@@ -261,7 +262,21 @@ double TFCS1DFunctionHistogram::linear(double y1,double y2,double x1,double x2,d
double n=y1-m*x1;
x=(y-n)/m;
}
+
+ return x;
+}
+double TFCS1DFunctionHistogram::non_linear(double y1,double y2,double x1,double x2,double y)
+{
+ double x=-1;
+ double eps=0.0000000001;
+ if((y2-y1)<eps) x=x1;
+ else
+ {
+ double b=(x1-x2)/(sqrt(y1)-sqrt(y2));
+ double a=x1-b*sqrt(y1);
+ x=a+b*sqrt(y);
+ }
return x;
}
@@ -271,39 +286,34 @@ double TFCS1DFunctionHistogram::get_inverse(double rnd)
double value = 0.;
- if(rnd<m_HistoContents[0]+(m_HistoContents[1]-m_HistoContents[0])/2.0)
+ if(rnd<m_HistoContents[0])
{
double x1=m_HistoBorders[0];
double x2=m_HistoBorders[1];
double y1=0;
- double y2=m_HistoContents[0]+(m_HistoContents[1]-m_HistoContents[0])/2.0;
- double x=linear(y1,y2,x1,x2,rnd);
+ double y2=m_HistoContents[0];
+ double x=non_linear(y1,y2,x1,x2,rnd);
value=x;
}
else
{
- for(unsigned int b=0;b<m_HistoContents.size()-1;b++)
+ //find the first HistoContent element that is larger than rnd:
+ vector<float>::iterator larger_element = std::upper_bound(m_HistoContents.begin(), m_HistoContents.end(), rnd);
+ int index=larger_element-m_HistoContents.begin();
+ double y=m_HistoContents[index];
+ double x1=m_HistoBorders[index];
+ double x2=m_HistoBorders[index+1];
+ double y1=m_HistoContents[index-1];
+ double y2=y;
+ if((index+1)==((int)m_HistoContents.size()-1))
{
- double y=m_HistoContents[b]+(m_HistoContents[b+1]-m_HistoContents[b])/2.0;
- if(y>rnd)
- {
- double x1=m_HistoBorders[b];
- double x2=m_HistoBorders[b+1];
- double y1=m_HistoContents[b-1]+(m_HistoContents[b]-m_HistoContents[b-1])/2.0;
- double y2=y;
- if((b+1)==m_HistoContents.size()-1)
- {
- x2=m_HistoBorders[m_HistoBorders.size()-1];
- y2=1;
- }
- double x=linear(y1,y2,x1,x2,rnd);
- value=x;
- break;
- }
+ x2=m_HistoBorders[m_HistoBorders.size()-1];
+ y2=m_HistoContents[m_HistoContents.size()-1];
}
+ double x=non_linear(y1,y2,x1,x2,rnd);
+ value=x;
}
return value;
}
-
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
index dc111f5ae897..60d9a85cf526 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
@@ -61,20 +61,15 @@ void TFCSPCAEnergyParametrization::simulate(TFCSSimulationState& simulstate,cons
{
int pcabin=simulstate.Ebin();
-
- TMatrixDSym* symCov =m_symCov[pcabin-1];
-
- TMatrixDSymEigen cov_eigen(*symCov);
- TMatrixD EV1 = cov_eigen.GetEigenVectors();
- TMatrixD* EV=&EV1;
-
- 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];
- TVectorD* LowerBounds =m_LowerBounds[pcabin-1];
+
+ 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];
+ 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;
@@ -87,7 +82,7 @@ void TFCSPCAEnergyParametrization::simulate(TFCSSimulationState& simulstate,cons
layer.push_back(thislayer);
}
layer.push_back("totalE");
-
+
double* vals_gauss_means=(double*)Gauss_means->GetMatrixArray();
double* vals_gauss_rms =Gauss_rms->GetMatrixArray();
double* vals_lowerBounds=LowerBounds->GetMatrixArray();
@@ -102,7 +97,7 @@ void TFCSPCAEnergyParametrization::simulate(TFCSSimulationState& simulstate,cons
double gauszz=Random3->Gaus(mean,rms);
input_data[l]=gauszz;
}
-
+
P2X(SigmaValues, MeanValues, EV, layer.size(), input_data, output_data, layer.size());
double *simdata_uniform = new double[layer.size()];
@@ -139,7 +134,8 @@ void TFCSPCAEnergyParametrization::simulate(TFCSSimulationState& simulstate,cons
}
}
- double total_energy=simdata[layer.size()-1];
+ double total_energy=simdata[layer.size()-1]*simulstate.E()/Ekin_nominal();
+ //double total_energy=simdata[layer.size()-1];
simulstate.set_E(total_energy);
ATH_MSG_DEBUG("set E to total_energy="<<total_energy);
@@ -181,13 +177,17 @@ void TFCSPCAEnergyParametrization::P2X(TVectorD* SigmaValues, TVectorD* MeanValu
}
}
}
-void TFCSPCAEnergyParametrization::loadInputs(TFile* file){
- loadInputs(file, "");
-}
-void TFCSPCAEnergyParametrization::loadInputs(TFile* file, std::string folder)
+bool TFCSPCAEnergyParametrization::loadInputs(TFile* file)
{
+ return loadInputs(file, "");
+}
+bool TFCSPCAEnergyParametrization::loadInputs(TFile* file, std::string folder)
+{
+
+ bool load_ok=1;
+
int trynext=1;
TString x;
if(folder=="") x="bin";
@@ -207,13 +207,17 @@ void TFCSPCAEnergyParametrization::loadInputs(TFile* file, std::string folder)
file->cd(x+"1/pca");
IntArray* RelevantLayers=(IntArray*)gDirectory->Get("RelevantLayers");
- if(RelevantLayers == NULL) {
+ if(RelevantLayers == NULL)
+ {
ATH_MSG_ERROR("TFCSPCAEnergyParametrization::m_RelevantLayers in first pcabin is null!");
- } else {
- m_RelevantLayers.reserve(RelevantLayers->GetSize());
- for(int i=0;i<RelevantLayers->GetSize();i++) m_RelevantLayers.push_back(RelevantLayers->GetAt(i));
+ load_ok=false;
}
-
+
+ if(!load_ok) return false;
+
+ m_RelevantLayers.reserve(RelevantLayers->GetSize());
+ for(int i=0;i<RelevantLayers->GetSize();i++) m_RelevantLayers.push_back(RelevantLayers->GetAt(i));
+
for(int bin=1;bin<=m_numberpcabins;bin++)
{
@@ -225,24 +229,28 @@ void TFCSPCAEnergyParametrization::loadInputs(TFile* file, std::string folder)
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::m_symCov in pcabin "<<bin<<" is null!");
- if(MeanValues == NULL) ATH_MSG_WARNING("TFCSPCAEnergyParametrization::m_MeanValues in pcabin "<<bin<<" is null!");
- if(SigmaValues == NULL) ATH_MSG_WARNING("TFCSPCAEnergyParametrization::m_SigmaValues in pcabin "<<bin<<" is null!");
- if(Gauss_means == NULL) ATH_MSG_WARNING("TFCSPCAEnergyParametrization::m_Gauss_means in pcabin "<<bin<<" is null!");
- if(Gauss_rms == NULL) ATH_MSG_WARNING("TFCSPCAEnergyParametrization::m_Gause_rms in pcabin "<<bin<<" is null!");
- if(LowerBounds == NULL) ATH_MSG_WARNING("TFCSPCAEnergyParametrization::m_LowerBounds in pcabin "<<bin<<" is null!");
-
- m_symCov.push_back(symCov);
+
+ 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;
+
+ TMatrixDSymEigen cov_eigen(*symCov);
+ TMatrixD *EV = new TMatrixD(cov_eigen.GetEigenVectors());
+ m_EV.push_back(EV);
m_MeanValues.push_back(MeanValues);
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;
-
+
for(unsigned int i=0;i<m_RelevantLayers.size();i++)
layerNr.push_back(m_RelevantLayers[i]);
@@ -273,4 +281,6 @@ void TFCSPCAEnergyParametrization::loadInputs(TFile* file, std::string folder)
}
+ return true;
+
}
--
GitLab