From 64a547771a4c6374f2671f6048d8678cf87ef889 Mon Sep 17 00:00:00 2001
From: Michael Duehrssen <michael.duehrssen@cern.ch>
Date: Thu, 26 Mar 2020 11:23:28 +0100
Subject: [PATCH] Address issue with negative hit energies seen in
 ATLASSIM-4474

---
 .../TFCSLateralShapeParametrizationFluctChain.cxx | 15 ++++++++++-----
 .../TFCSLateralShapeParametrizationHitChain.cxx   | 12 ++++++++----
 2 files changed, 18 insertions(+), 9 deletions(-)

diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationFluctChain.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationFluctChain.cxx
index 26b99bbab8f..a45ea62759d 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationFluctChain.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationFluctChain.cxx
@@ -35,6 +35,12 @@ float TFCSLateralShapeParametrizationFluctChain::get_E_hit(TFCSSimulationState&
 
 FCSReturnCode TFCSLateralShapeParametrizationFluctChain::simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const
 {
+  const float Elayer=simulstate.E(calosample());
+  if (Elayer == 0) {
+    ATH_MSG_VERBOSE("TFCSLateralShapeParametrizationHitChain::simulate(): Elayer=0, nothing to do");
+    return FCSSuccess;
+  }
+  
   // Call get_number_of_hits() only once, as it could contain a random number
   float sigma2  = get_sigma2_fluctuation(simulstate, truth, extrapol);
   if (sigma2 >= s_max_sigma2_fluctuation) {
@@ -45,11 +51,9 @@ FCSReturnCode TFCSLateralShapeParametrizationFluctChain::simulate(TFCSSimulation
   //Limit to relative precision of 10^-4=0.01%. ATLAS calorimeter is ~0.1% at best
   if(sigma2<1e-8) sigma2=1e-8; 
 
-  const float Elayer=simulstate.E(calosample());
-
   //Make a good guess of the needed hit energy, assuming all hits would have the same energy
   const float Eavghit=get_E_hit(simulstate,truth,extrapol);
-  const float Eavghit_tenth=Eavghit/10;
+  const float absEavghit_tenth=fabs(Eavghit/10);
   float sumEhit=0;
   float error2_sumEhit=0;
   float error2=2*s_max_sigma2_fluctuation;
@@ -69,7 +73,7 @@ FCSReturnCode TFCSLateralShapeParametrizationFluctChain::simulate(TFCSSimulation
     //hit.E()=Eavghit;
     do {
       hit.E()=CLHEP::RandGauss::shoot(simulstate.randomEngine(), Eavghit, m_RMS*Eavghit);
-    } while (hit.E()<Eavghit_tenth);
+    } while (fabs(hit.E())<absEavghit_tenth);
     bool failed=false;
     for(TFCSLateralShapeParametrizationHitBase* hitsim : m_chain) {
       if (debug) {
@@ -91,8 +95,9 @@ FCSReturnCode TFCSLateralShapeParametrizationFluctChain::simulate(TFCSSimulation
       ++ihit;
       const float Ehit=hit.E();
       sumEhit+=Ehit;
+      float sumEhit2=sumEhit*sumEhit;
       error2_sumEhit+=Ehit*Ehit;
-      if(sumEhit>0) error2=error2_sumEhit/(sumEhit*sumEhit);
+      if(sumEhit2>0) error2=error2_sumEhit/sumEhit2;
     } else {
       if (ifail >= FCS_RETRY_COUNT) {
         ATH_MSG_ERROR("TFCSLateralShapeParametrizationFluctChain::simulate(): simulate_hit call failed after " << FCS_RETRY_COUNT << "retries");
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitChain.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitChain.cxx
index 004a900104f..a7c615c83b1 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitChain.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitChain.cxx
@@ -91,9 +91,13 @@ float TFCSLateralShapeParametrizationHitChain::get_sigma2_fluctuation(TFCSSimula
 FCSReturnCode TFCSLateralShapeParametrizationHitChain::simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const
 {
   const float Elayer=simulstate.E(calosample());
+  if (Elayer == 0) {
+    ATH_MSG_VERBOSE("TFCSLateralShapeParametrizationHitChain::simulate(): Elayer=0, nothing to do");
+    return FCSSuccess;
+  }
   const float Ehit=get_E_hit(simulstate,truth,extrapol);
-  if (Ehit <= 0) {
-    ATH_MSG_ERROR("TFCSLateralShapeParametrizationHitChain::simulate(): Ehit negative Ehit="<<Ehit);
+  if (Ehit * Elayer <= 0) {
+    ATH_MSG_ERROR("TFCSLateralShapeParametrizationHitChain::simulate(): Ehit and Elayer have different sign! Ehit="<<Ehit<<" Elayer="<<Elayer);
     return FCSFatal;
   }
 
@@ -151,13 +155,13 @@ FCSReturnCode TFCSLateralShapeParametrizationHitChain::simulate(TFCSSimulationSt
     ++ihit;
     
     if( (ihit==10*nhit) || (ihit==100*nhit) ) {
-      ATH_MSG_WARNING("TFCSLateralShapeParametrizationHitChain::simulate(): Iterated " << ihit << " times, expected " << nhit <<" times. Deposited E("<<calosample()<<")="<<sumEhit);
+      ATH_MSG_WARNING("TFCSLateralShapeParametrizationHitChain::simulate(): Iterated " << ihit << " times, expected " << nhit <<" times. Deposited E("<<calosample()<<")="<<sumEhit<<" expected E="<<Elayer);
     }                                                                                                                         
     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);
+  } while (fabs(sumEhit)<fabs(Elayer));
 
   if (debug) {
     for(TFCSLateralShapeParametrizationHitBase* reset : m_chain) {
-- 
GitLab