From 9370e97a5d361670efa899dfaab2f0999f199e06 Mon Sep 17 00:00:00 2001
From: Petr Jacka <petr.jacka@cern.ch>
Date: Wed, 8 Jan 2020 15:03:13 +0000
Subject: [PATCH] Using get_E_hit function in
TFCSLateralShapeParametrizationFluctChain and
TFCSLateralShapeParametrizationHitChain
---
.../TFCSHistoLateralShapeWeight.h | 9 ++++-
...FCSLateralShapeParametrizationFluctChain.h | 3 ++
.../TFCSLateralShapeParametrizationHitBase.h | 7 ++++
.../TFCSLateralShapeParametrizationHitChain.h | 7 ++++
.../src/TFCSHistoLateralShapeWeight.cxx | 15 +++++--
.../TFCSHistoLateralShapeWeightHitAndMiss.cxx | 37 ++++++++++++-----
...SLateralShapeParametrizationFluctChain.cxx | 15 ++++++-
...TFCSLateralShapeParametrizationHitBase.cxx | 18 +++++++++
...FCSLateralShapeParametrizationHitChain.cxx | 40 +++++++++++++++++--
9 files changed, 132 insertions(+), 19 deletions(-)
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h
index dc650c2ca08d4..3af555f21ce84 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h
@@ -21,12 +21,19 @@ public:
bool Initialize(TH1* hist);
virtual void Print(Option_t *option = "") const override;
+ virtual void setMinWeight(float minWeight) {m_minWeight=minWeight;}
+ virtual void setMaxWeight(float maxWeight) {m_maxWeight=maxWeight;}
+ virtual float getMinWeight() const override;
+ virtual float getMaxWeight() const override;
+
protected:
/// Histogram to be used for the shape simulation
/// The histogram x-axis should be in dR^2=deta^2+dphi^2
TH1* m_hist{nullptr};
+ float m_minWeight{-1.};
+ float m_maxWeight{-1.};
- ClassDefOverride(TFCSHistoLateralShapeWeight,1) //TFCSHistoLateralShapeWeight
+ ClassDefOverride(TFCSHistoLateralShapeWeight,2) //TFCSHistoLateralShapeWeight
};
#if defined(__ROOTCLING__) && defined(__FastCaloSimStandAlone__)
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationFluctChain.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationFluctChain.h
index c8e218c3c19e1..1d94eb14e8dee 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationFluctChain.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationFluctChain.h
@@ -15,6 +15,9 @@ public:
// set and get the amount of hit energy fluctation around E/n for n hits
void set_hit_RMS(float RMS);
float get_hit_RMS() const {return m_RMS;};
+
+ ///Get hit energy from layer energy and number of hits
+ virtual float get_E_hit(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const override;
virtual FCSReturnCode simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const override;
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitBase.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitBase.h
index 44635409f11c0..f030ca57243bc 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitBase.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitBase.h
@@ -18,6 +18,13 @@ public:
///Call get_number_of_hits() only once per shower simulation, as it could be calculated with random numbers and give different results each time. Return a value of -1 if this instance can't determine
virtual int get_number_of_hits(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const;
+
+ ///Get hit energy from layer energy and number of hits
+ virtual float get_E_hit(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const;
+
+ ///Get minimum and maximum value of weight for hit energy reweighting
+ virtual float getMinWeight() const;
+ virtual float getMaxWeight() const;
class Hit
{
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitChain.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitChain.h
index 57027b4604019..0717b2ed338fb 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitChain.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitChain.h
@@ -33,9 +33,16 @@ public:
/// Call get_number_of_hits() only once, as it could contain a random number
virtual int get_number_of_hits(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const;
+
+ ///Get hit energy from layer energy and number of hits
+ virtual float get_E_hit(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const;
///Give the effective size sigma^2 of the fluctuations that should be generated.
virtual float get_sigma2_fluctuation(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const;
+
+ ///Get minimum and maximum value of weight for hit energy reweighting
+ virtual float getMinWeight() const;
+ virtual float getMaxWeight() const;
static constexpr float s_max_sigma2_fluctuation=1000;//! Do not persistify!
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeWeight.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeWeight.cxx
index 8cee173ab4721..9eee5d468cef4 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeWeight.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeWeight.cxx
@@ -18,14 +18,23 @@
TFCSHistoLateralShapeWeight::TFCSHistoLateralShapeWeight(const char* name, const char* title) :
TFCSLateralShapeParametrizationHitBase(name,title)
-{
-}
+{}
TFCSHistoLateralShapeWeight::~TFCSHistoLateralShapeWeight()
{
if(m_hist) delete m_hist;
}
+float TFCSHistoLateralShapeWeight::getMinWeight() const
+{
+ return m_minWeight;
+}
+
+float TFCSHistoLateralShapeWeight::getMaxWeight() const
+{
+ return m_maxWeight;
+}
+
FCSReturnCode TFCSHistoLateralShapeWeight::simulate_hit(Hit& hit,TFCSSimulationState& simulstate,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* /*extrapol*/)
{
if (!simulstate.randomEngine()) {
@@ -82,7 +91,7 @@ void TFCSHistoLateralShapeWeight::Print(Option_t *option) const
TFCSLateralShapeParametrizationHitBase::Print(option);
if(longprint) {
- if(m_hist) ATH_MSG_INFO(optprint <<" Histogram: "<<m_hist->GetNbinsX()<<" bins ["<<m_hist->GetXaxis()->GetXmin()<<","<<m_hist->GetXaxis()->GetXmax()<<"]");
+ if(m_hist) ATH_MSG_INFO(optprint <<" Histogram: "<<m_hist->GetNbinsX()<<" bins ["<<m_hist->GetXaxis()->GetXmin()<<","<<m_hist->GetXaxis()->GetXmax()<<"]" << " min weight: " << m_minWeight << " max weight: " << m_maxWeight);
else ATH_MSG_INFO(optprint <<" no Histogram");
}
}
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeWeightHitAndMiss.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeWeightHitAndMiss.cxx
index 735529227470c..608d1731ed040 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeWeightHitAndMiss.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeWeightHitAndMiss.cxx
@@ -55,20 +55,37 @@ FCSReturnCode TFCSHistoLateralShapeWeightHitAndMiss::simulate_hit(Hit& hit,TFCSS
if(RMS>0) {
weight=CLHEP::RandGauss::shoot(simulstate.randomEngine(), meanweight, RMS);
}
- if(meanweight<=1) {
- //if meanweight<=1, give lower energy to hit.
- //TFCSLateralShapeParametrizationHitChain needs to be able to generate more hits in this case
- hit.E()*=weight;
- } else {
- //if meanweight>1, accept only 1/meanweight events, but give them a higher energy increased by weight.
- //This leads to larger fluctuations, while keeping the shape unchanged.
- float prob=1.0/meanweight;
+
+ /* -------------------------------------------------------------------
+ * Weight is used to scale hit energy.
+ *
+ * if (meanweight > m_minWeight and meanweight < m_maxWeight)
+ * Hit is accecpted with probability of m_minWeight/meanweight.
+ * If not accepted, weight is set to zero (this is
+ * equivalent to not accept the hit).
+ *
+ * if (meanweight<=m_minWeight)
+ * Give lower energy to hit with probability 1.
+ * TFCSLateralShapeParametrizationHitChain needs to be able
+ * to generate more hits in this case
+ *
+ * if (meanweight >= m_maxWeight)
+ * meanweight is above upper threshold. It is set to m_maxWeight.
+ * -------------------------------------------------------------------
+ */
+ if (meanweight >= m_maxWeight) {
+ weight = m_maxWeight;
+ }
+ else if (meanweight > m_minWeight){
+ float prob=m_minWeight/meanweight;
float rnd=CLHEP::RandFlat::shoot(simulstate.randomEngine());
- if(rnd<prob) hit.E()*=weight;
- else hit.E()=0;
+ if(rnd>=prob) weight = 0.;
}
+
+ hit.E()*=weight;
ATH_MSG_DEBUG("HIT: E="<<hit.E()<<" dR_mm="<<delta_r_mm<<" meanweight="<<meanweight<<" weight="<<weight);
+
return FCSSuccess;
}
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationFluctChain.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationFluctChain.cxx
index 7f4eb1c71e453..26b99bbab8fa5 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationFluctChain.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationFluctChain.cxx
@@ -21,6 +21,18 @@ TFCSLateralShapeParametrizationFluctChain::TFCSLateralShapeParametrizationFluctC
{
}
+float TFCSLateralShapeParametrizationFluctChain::get_E_hit(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const
+{
+ const float sigma2 = get_sigma2_fluctuation(simulstate,truth,extrapol);
+ const int sample = calosample();
+ if(sigma2<=0 || sample<0) return -1.;
+ const float maxWeight = getMaxWeight();
+
+ if(maxWeight>0) return simulstate.E(sample)*sigma2/maxWeight;
+ else return simulstate.E(sample)*sigma2;
+
+}
+
FCSReturnCode TFCSLateralShapeParametrizationFluctChain::simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const
{
// Call get_number_of_hits() only once, as it could contain a random number
@@ -34,10 +46,9 @@ FCSReturnCode TFCSLateralShapeParametrizationFluctChain::simulate(TFCSSimulation
if(sigma2<1e-8) sigma2=1e-8;
const float Elayer=simulstate.E(calosample());
- const int nhit = 1.0 / sigma2;
//Make a good guess of the needed hit energy, assuming all hits would have the same energy
- const float Eavghit=Elayer/nhit;
+ const float Eavghit=get_E_hit(simulstate,truth,extrapol);
const float Eavghit_tenth=Eavghit/10;
float sumEhit=0;
float error2_sumEhit=0;
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitBase.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitBase.cxx
index 8e2badda702e5..e2b2f33af8239 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitBase.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitBase.cxx
@@ -26,6 +26,24 @@ int TFCSLateralShapeParametrizationHitBase::get_number_of_hits(TFCSSimulationSta
return -1;
}
+float TFCSLateralShapeParametrizationHitBase::get_E_hit(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const
+{
+ const int nhits = get_number_of_hits(simulstate,truth,extrapol);
+ const int sample = calosample();
+ if(nhits<=0 || sample<0) return -1.;
+ else return simulstate.E(sample)/nhits;
+}
+
+float TFCSLateralShapeParametrizationHitBase::getMinWeight() const
+{
+ return -1.;
+}
+
+float TFCSLateralShapeParametrizationHitBase::getMaxWeight() const
+{
+ return -1.;
+}
+
FCSReturnCode TFCSLateralShapeParametrizationHitBase::simulate_hit(Hit& hit,TFCSSimulationState& /*simulstate*/,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* extrapol)
{
int cs=calosample();
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitChain.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitChain.cxx
index e1d27921966a7..c0fe229ae2dc3 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitChain.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitChain.cxx
@@ -43,6 +43,36 @@ int TFCSLateralShapeParametrizationHitChain::get_number_of_hits(TFCSSimulationSt
return 1;
}
+float TFCSLateralShapeParametrizationHitChain::getMinWeight() const
+{
+ for(TFCSLateralShapeParametrizationHitBase* hitsim : m_chain) {
+ float weight=hitsim->getMinWeight();
+ if(weight>0.) return weight;
+ }
+ return -1.;
+}
+
+float TFCSLateralShapeParametrizationHitChain::getMaxWeight() const
+{
+ for(TFCSLateralShapeParametrizationHitBase* hitsim : m_chain) {
+ float weight=hitsim->getMaxWeight();
+ if(weight>0.) return weight;
+ }
+ return -1.;
+}
+
+
+float TFCSLateralShapeParametrizationHitChain::get_E_hit(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const
+{
+ const int nhits = get_number_of_hits(simulstate,truth,extrapol);
+ const int sample = calosample();
+ if(nhits<=0 || sample<0) return -1.;
+ const float maxWeight = getMaxWeight();// maxWeight = -1 if shapeWeight class is not in m_chain
+
+ if(maxWeight>0) return simulstate.E(sample)/(maxWeight*nhits); // maxWeight is used only if shapeWeight class is in m_chain
+ else return simulstate.E(sample)/nhits; // Otherwise, old definition of E_hit is used
+}
+
float TFCSLateralShapeParametrizationHitChain::get_sigma2_fluctuation(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const
{
if(m_number_of_hits_simul) {
@@ -67,7 +97,7 @@ FCSReturnCode TFCSLateralShapeParametrizationHitChain::simulate(TFCSSimulationSt
}
const float Elayer=simulstate.E(calosample());
- const float Ehit=Elayer/nhit;
+ const float Ehit=get_E_hit(simulstate,truth,extrapol);
float sumEhit=0;
const bool debug = msgLvl(MSG::DEBUG);
@@ -107,8 +137,12 @@ FCSReturnCode TFCSLateralShapeParametrizationHitChain::simulate(TFCSSimulationSt
}
sumEhit+=hit.E();
++ihit;
- if(ihit>100*nhit && ihit>1000) {
- ATH_MSG_WARNING("TFCSLateralShapeParametrizationHitChain::simulate(): aborting hit chain, iterated " << 100*nhit << " times, expected " << nhit<<" times. Deposited E("<<calosample()<<")="<<sumEhit<<" expected E="<<Elayer);
+
+ if( (ihit==10*nhit) || (ihit==100*nhit) ) {
+ ATH_MSG_WARNING("TFCSLateralShapeParametrizationHitChain::simulate(): Iterated " << ihit << " times, expected " << nhit <<" times. Deposited E("<<calosample()<<")="<<sumEhit);
+ }
+ if(ihit>1000*nhit && ihit>1000) {
+ ATH_MSG_WARNING("TFCSLateralShapeParametrizationHitChain::simulate(): Aborting hit chain, iterated " << 1000*nhit << " times, expected " << nhit <<" times. Deposited E("<<calosample()<<")="<<sumEhit<<" expected E="<<Elayer);
break;
}
} while (sumEhit<Elayer);
--
GitLab