diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSHitCellMappingWiggleEMB.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSHitCellMappingWiggleEMB.h
index 1624d6a19d6fa6352c7301de8f18719aa413b02e..e776860fd7275ae7f6e8c45250674b2582a70af1 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSHitCellMappingWiggleEMB.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSHitCellMappingWiggleEMB.h
@@ -21,7 +21,7 @@ private:
   double m_wiggleLayer2[50];
   double m_wiggleLayer3[50];
 
-  double doWiggle();
+  double doWiggle(double searchRand);
 
   ClassDefOverride(TFCSHitCellMappingWiggleEMB,1)  //TFCSHitCellMappingWiggleEMB
 };
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitNumberFromE.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitNumberFromE.h
index ec4ef2f20d365d47f7485adadb87b4c8014ff7e1..cbe0e39996096567ab2a9130c8647e5eb5a545bb 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitNumberFromE.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitNumberFromE.h
@@ -8,7 +8,6 @@
 #include "ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitBase.h"
 
 #include "TH2.h"
-#include "TRandom3.h"
 
 
 class TFCSLateralShapeParametrizationHitNumberFromE:public TFCSLateralShapeParametrizationHitBase {
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 f2ef1a4dc74ad6c7c71f93003bd2b8a6ee8204e4..a8b8fe39f9e9390fad9b92e582a3d79cfce542c6 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSSimulationState.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSSimulationState.h
@@ -12,10 +12,18 @@
 
 class CaloDetDescrElement;
 
+namespace CLHEP
+{
+ class HepRandomEngine;
+}
+
 class TFCSSimulationState:public TObject
 {
   public:
-    TFCSSimulationState();
+    TFCSSimulationState(CLHEP::HepRandomEngine* randomEngine = nullptr);
+
+    CLHEP::HepRandomEngine* randomEngine() { return m_randomEngine; }
+    void   setRandomEngine(CLHEP::HepRandomEngine *engine) { m_randomEngine = engine; }
 
     bool   is_valid() const {return m_Ebin>=0;};
     double E() const {return m_Etot;};
@@ -41,6 +49,8 @@ class TFCSSimulationState:public TObject
 
     void clear();
   private:
+    CLHEP::HepRandomEngine* m_randomEngine;
+      
     int    m_Ebin;
     double m_Etot;
     // TO BE CLEANED UP! SHOULD ONLY STORE EITHER E OR EFRAC!!!
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyBinParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyBinParametrization.cxx
index ba3db006ec7baa9a30349d5ce2d1b5adcbada951..9938530a397817cd1a7f98def9447640388ed9ce 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyBinParametrization.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSEnergyBinParametrization.cxx
@@ -2,11 +2,13 @@
   Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
 */
 
+#include "CLHEP/Random/RandFlat.h"
+
 #include "ISF_FastCaloSimEvent/TFCSEnergyBinParametrization.h"
 #include "ISF_FastCaloSimEvent/TFCSTruthState.h"
 #include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
+
 #include "TMath.h"
-#include "TRandom.h"
 
 //=============================================
 //======= TFCSSelectEnergyBin =========
@@ -101,12 +103,16 @@ void TFCSEnergyBinParametrization::Print(Option_t *option) const
 
 FCSReturnCode TFCSEnergyBinParametrization::simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* /*extrapol*/)
 {
+  if (!simulstate.randomEngine()) {
+    return FCSFatal;
+  }
+
   int pdgid=truth->pdgid();
   if(!is_match_pdgid(pdgid)) {
     ATH_MSG_ERROR("TFCSEnergyBinParametrization::simulate(): cannot simulate pdgid="<<pdgid);
     return FCSFatal;
   }
-  float searchRand=gRandom->Rndm();
+  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()) {
     ATH_MSG_ERROR("TFCSEnergyBinParametrization::simulate(): cannot simulate bin="<<chosenBin);
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeParametrization.cxx
index d726ffd45729fcb4c5ad8f836113011127b71036..fe310de94c0f6a2dd1948b9efdf2c64c77c8cdbf 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeParametrization.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeParametrization.cxx
@@ -2,6 +2,9 @@
   Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
 */
 
+#include "CLHEP/Random/RandFlat.h"
+#include "CLHEP/Random/RandPoisson.h"
+
 #include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeParametrization.h"
 #include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h"
 #include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
@@ -9,7 +12,6 @@
 
 #include "TFile.h"
 #include "TMath.h"
-#include "TRandom3.h"
 #include "TH2.h"
 
 
@@ -27,9 +29,13 @@ TFCSHistoLateralShapeParametrization::~TFCSHistoLateralShapeParametrization()
 {
 }
 
-int TFCSHistoLateralShapeParametrization::get_number_of_hits(TFCSSimulationState& /*simulstate*/,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* /*extrapol*/) const
+int TFCSHistoLateralShapeParametrization::get_number_of_hits(TFCSSimulationState &simulstate, const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* /*extrapol*/) const
 {
-  return gRandom->Poisson(m_nhits);
+  if (!simulstate.randomEngine()) {
+    return -1;
+  }
+
+  return CLHEP::RandPoisson::shoot(simulstate.randomEngine(), m_nhits);
 }
 
 void TFCSHistoLateralShapeParametrization::set_number_of_hits(float nhits)
@@ -37,8 +43,12 @@ void TFCSHistoLateralShapeParametrization::set_number_of_hits(float nhits)
   m_nhits=nhits;
 }
 
-FCSReturnCode TFCSHistoLateralShapeParametrization::simulate_hit(Hit& hit,TFCSSimulationState& /*simulstate*/,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* extrapol)
+FCSReturnCode TFCSHistoLateralShapeParametrization::simulate_hit(Hit &hit, TFCSSimulationState &simulstate, const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* extrapol)
 {
+  if (!simulstate.randomEngine()) {
+    return FCSFatal;
+  }
+
   const int cs=calosample();
   const double center_eta=0.5*( extrapol->eta(cs, CaloSubPos::SUBPOS_ENT) + extrapol->eta(cs, CaloSubPos::SUBPOS_EXT) );
   const double center_phi=0.5*( extrapol->phi(cs, CaloSubPos::SUBPOS_ENT) + extrapol->phi(cs, CaloSubPos::SUBPOS_EXT) );
@@ -46,10 +56,8 @@ FCSReturnCode TFCSHistoLateralShapeParametrization::simulate_hit(Hit& hit,TFCSSi
   const double center_z=0.5*( extrapol->z(cs, CaloSubPos::SUBPOS_ENT) + extrapol->z(cs, CaloSubPos::SUBPOS_EXT) );
 
   float alpha, r, rnd1, rnd2;
-  //The use of 1-gRandom->Rndm() is a fudge for TRandom3, as it gererates random numbers in (0,1], but [0,1) or (0,1) is needed. 
-  //CLHEP should generate random numbers in (0,1), so this fudge is no longer needed after migrating to CLHEP random numbers
-  rnd1=1-gRandom->Rndm(); 
-  rnd2=1-gRandom->Rndm();
+  rnd1 = CLHEP::RandFlat::shoot(simulstate.randomEngine());
+  rnd2 = CLHEP::RandFlat::shoot(simulstate.randomEngine());
   if(is_phi_symmetric()) {
     if(rnd2>=0.5) { //Fill negative phi half of shape
       rnd2-=0.5;
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHitCellMappingWiggle.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHitCellMappingWiggle.cxx
index 70f8cf2f0e779933b8bde30b34f08cf939790236..b08afd6711c64923e059c79790900b48fc07cd03 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHitCellMappingWiggle.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHitCellMappingWiggle.cxx
@@ -2,13 +2,15 @@
   Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
 */
 
+#include "CLHEP/Random/RandFlat.h"
+
 #include "ISF_FastCaloSimEvent/TFCSHitCellMappingWiggle.h"
 #include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
 #include "ISF_FastCaloSimEvent/TFCSTruthState.h"
 #include "ISF_FastCaloSimEvent/TFCSExtrapolationState.h"
 #include "ISF_FastCaloSimEvent/TFCS1DFunctionInt32Histogram.h"
+
 #include "TH1.h"
-#include "TRandom3.h"
 #include "TVector2.h"
 #include "TMath.h"
 
@@ -87,6 +89,10 @@ void TFCSHitCellMappingWiggle::initialize(const std::vector< const TH1* > histog
 
 FCSReturnCode TFCSHitCellMappingWiggle::simulate_hit(Hit& hit,TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol)
 {
+  if (!simulstate.randomEngine()) {
+    return FCSFatal;
+  }
+
   float eta=fabs(hit.eta());
   if(eta<m_bin_low_edge[0] || eta>=m_bin_low_edge[get_number_of_bins()]) {
     return TFCSHitCellMapping::simulate_hit(hit,simulstate,truth,extrapol);
@@ -97,8 +103,8 @@ FCSReturnCode TFCSHitCellMappingWiggle::simulate_hit(Hit& hit,TFCSSimulationStat
 
   const TFCS1DFunction* func=get_function(bin);
   if(func) {
-    double rnd = 1-gRandom->Rndm();
-    
+    double rnd = CLHEP::RandFlat::shoot(simulstate.randomEngine());
+
     double wiggle=func->rnd_to_fct(rnd);
 
     ATH_MSG_DEBUG("HIT: E="<<hit.E()<<" cs="<<calosample()<<" eta="<<hit.eta()<<" phi="<<hit.phi()<<" wiggle="<<wiggle<<" bin="<<bin<<" ["<<get_bin_low_edge(bin)<<","<<get_bin_up_edge(bin)<<"] func="<<func);
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHitCellMappingWiggleEMB.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHitCellMappingWiggleEMB.cxx
index b817853ee81ad3e37eaf1b0d4bc4570f5e7e7a1b..912250b89d950fc0f100ce75d5ae459971524af0 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHitCellMappingWiggleEMB.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHitCellMappingWiggleEMB.cxx
@@ -2,9 +2,11 @@
   Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
 */
 
+#include "CLHEP/Random/RandFlat.h"
+
 #include "ISF_FastCaloSimEvent/TFCSHitCellMappingWiggleEMB.h"
 #include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
-#include "TRandom3.h"
+
 #include "TVector2.h"
 #include "TMath.h"
 
@@ -27,7 +29,7 @@ TFCSHitCellMappingWiggleEMB::TFCSHitCellMappingWiggleEMB(const char* name, const
   }
 }
 
-double TFCSHitCellMappingWiggleEMB::doWiggle()
+double TFCSHitCellMappingWiggleEMB::doWiggle(double searchRand)
 {
  int layer=calosample();
  
@@ -44,10 +46,7 @@ double TFCSHitCellMappingWiggleEMB::doWiggle()
  {
   return 0.0;
  }
- 
- //Set random numbers
- double searchRand = gRandom->Rndm();
-  
+
  //Now for layer dependant approach
  if(layer == 1)
  {
@@ -75,10 +74,14 @@ double TFCSHitCellMappingWiggleEMB::doWiggle()
 
 FCSReturnCode TFCSHitCellMappingWiggleEMB::simulate_hit(Hit& hit,TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol)
 {
+  if (!simulstate.randomEngine()) {
+    return FCSFatal;
+  }
+
   int cs=calosample();
 
   double wiggle = 0.0;
-  if(cs < 4 && cs > 0) wiggle = doWiggle();
+  if(cs < 4 && cs > 0) wiggle = doWiggle(CLHEP::RandFlat::shoot(simulstate.randomEngine()));
 
   ATH_MSG_DEBUG("HIT: E="<<hit.E()<<" cs="<<cs<<" eta="<<hit.eta()<<" phi="<<hit.phi()<<" wiggle="<<wiggle);
 
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitChain.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitChain.cxx
index 798689be0fd32219b3700e00c8dccc9cde50da9f..7e3e400a63bb1a052c92b8d8e7956a56a40c0894 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitChain.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitChain.cxx
@@ -30,6 +30,7 @@ void TFCSLateralShapeParametrizationHitChain::set_geometry(ICaloGeometry* geo)
 
 int TFCSLateralShapeParametrizationHitChain::get_number_of_hits(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const
 {
+  // TODO: should we still do it?
   if(m_number_of_hits_simul) {
     int n=m_number_of_hits_simul->get_number_of_hits(simulstate,truth,extrapol);
     if(n<1) n=1;
@@ -45,7 +46,12 @@ int TFCSLateralShapeParametrizationHitChain::get_number_of_hits(TFCSSimulationSt
 FCSReturnCode TFCSLateralShapeParametrizationHitChain::simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol)
 {
   // Call get_number_of_hits() only once, as it could contain a random number
-  int nhit=get_number_of_hits(simulstate,truth,extrapol);
+  int nhit = get_number_of_hits(simulstate, truth, extrapol);
+  if (nhit <= 0) {
+    ATH_MSG_ERROR("TFCSLateralShapeParametrizationHitChain::simulate(): number of hits could not be calculated");
+    return FCSFatal;
+  }
+
   float Ehit=simulstate.E(calosample())/nhit;
 
   bool debug = msgLvl(MSG::DEBUG);
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitNumberFromE.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitNumberFromE.cxx
index ed0f9047fbc4e06946d7558db246bcdcc6c48eff..cbc4b4f2ad7f493c48d42efa027bac8679beda2d 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitNumberFromE.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitNumberFromE.cxx
@@ -2,6 +2,8 @@
   Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
 */
 
+#include "CLHEP/Random/RandPoisson.h"
+
 #include "ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitNumberFromE.h"
 #include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
 
@@ -19,11 +21,15 @@ TFCSLateralShapeParametrizationHitNumberFromE::TFCSLateralShapeParametrizationHi
 
 int TFCSLateralShapeParametrizationHitNumberFromE::get_number_of_hits(TFCSSimulationState& simulstate,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* /*extrapol*/) const
 {
+  if (!simulstate.randomEngine()) {
+    return -1;
+  }
+
   int cs=calosample();
   double energy=simulstate.E(cs);
 
   double sigma_stochastic=m_stochastic/sqrt(energy/1000.0);
-  int hits = gRandom->Poisson(1.0 / (sigma_stochastic*sigma_stochastic + m_constant*m_constant));
+  int hits = CLHEP::RandPoisson::shoot(simulstate.randomEngine(), 1.0 / (sigma_stochastic*sigma_stochastic + m_constant*m_constant));
 
   ATH_MSG_DEBUG("#hits="<<hits);
   
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
index 29a2e44b841f7b65bd9d93f9ce65f9665f38ae10..cfe5a29927dbab5247231500107f698ad4a9d64e 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
@@ -2,6 +2,8 @@
   Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
 */
 
+#include "CLHEP/Random/RandGauss.h"
+
 #include "ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h"
 #include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h"
 
@@ -11,7 +13,6 @@
 #include "TFile.h"
 #include "TKey.h"
 #include "TClass.h"
-#include "TRandom3.h"
 #include "TMatrixD.h"
 #include "TMatrixDSymEigen.h"
 #include "TMath.h"
@@ -60,6 +61,9 @@ 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();
  
@@ -69,8 +73,6 @@ FCSReturnCode TFCSPCAEnergyParametrization::simulate(TFCSSimulationState& simuls
   TVectorD* Gauss_means   =m_Gauss_means[pcabin-1];
   TVectorD* Gauss_rms     =m_Gauss_rms[pcabin-1];
   std::vector<TFCS1DFunction*> cumulative=m_cumulative[pcabin-1];
-  
-  TRandom3* Random3=new TRandom3(); Random3->SetSeed(0);
 
   std::vector<int> layerNr;
   for(unsigned int i=0;i<m_RelevantLayers.size();i++)
@@ -86,7 +88,7 @@ FCSReturnCode TFCSPCAEnergyParametrization::simulate(TFCSSimulationState& simuls
   {
    double mean=vals_gauss_means[l];
    double rms =vals_gauss_rms[l];
-   double gauszz=Random3->Gaus(mean,rms);
+   double gauszz = CLHEP::RandGauss::shoot(simulstate.randomEngine(), mean, rms);
    input_data[l]=gauszz;
   }
 
@@ -129,7 +131,6 @@ FCSReturnCode TFCSPCAEnergyParametrization::simulate(TFCSSimulationState& simuls
    simulstate.set_SF(scalefactor);
   }
 
-  delete Random3;
   delete [] output_data;
   delete [] input_data;
   delete [] simdata;
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrizationEkinSelectChain.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrizationEkinSelectChain.cxx
index cd01c674d62c0af960051315d9ef999cf2b0c7ce..9ecf2a23e8b9d01c607e9749e38d5866bac58b23 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrizationEkinSelectChain.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrizationEkinSelectChain.cxx
@@ -2,6 +2,8 @@
   Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
 */
 
+#include "CLHEP/Random/RandFlat.h"
+
 #include "ISF_FastCaloSimEvent/TFCSParametrizationEkinSelectChain.h"
 #include "ISF_FastCaloSimEvent/TFCSInvisibleParametrization.h"
 #include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
@@ -9,7 +11,6 @@
 #include "ISF_FastCaloSimEvent/TFCSExtrapolationState.h"
 
 #include <iostream>
-#include "TRandom.h"
 
 //=============================================
 //======= TFCSParametrizationEkinSelectChain =========
@@ -32,8 +33,9 @@ void TFCSParametrizationEkinSelectChain::push_back_in_bin(TFCSParametrizationBas
   push_back_in_bin(param,param->Ekin_min(),param->Ekin_max());
 }
 
-int TFCSParametrizationEkinSelectChain::get_bin(TFCSSimulationState&,const TFCSTruthState* truth, const TFCSExtrapolationState*) const
+int TFCSParametrizationEkinSelectChain::get_bin(TFCSSimulationState &simulstate, const TFCSTruthState* truth, const TFCSExtrapolationState*) const
 {
+  // TODO: fail for random engines
   float Ekin=truth->Ekin();
   int bin=val_to_bin(Ekin);
   
@@ -64,7 +66,7 @@ int TFCSParametrizationEkinSelectChain::get_bin(TFCSSimulationState&,const TFCST
     float denominator=logEkin_nominal-logEkin_previous;
     if(denominator<=0) return bin;
 
-    float rnd=gRandom->Rndm();
+    float rnd = CLHEP::RandFlat::shoot(simulstate.randomEngine());
     if(numerator/denominator<rnd) bin=prevbin;
     ATH_MSG_DEBUG("logEkin="<<logEkin<<" logEkin_previous="<<logEkin_previous<<" logEkin_nominal="<<logEkin_nominal<<" (rnd="<<1-rnd<<" < p(previous)="<<(1-numerator/denominator)<<")? => orgbin="<<prevbin+1<<" selbin="<<bin);
   } else {
@@ -83,7 +85,7 @@ int TFCSParametrizationEkinSelectChain::get_bin(TFCSSimulationState&,const TFCST
     float denominator=logEkin_next-logEkin_nominal;
     if(denominator<=0) return bin;
 
-    float rnd=gRandom->Rndm();
+    float rnd = CLHEP::RandFlat::shoot(simulstate.randomEngine());
     if(rnd<numerator/denominator) bin=nextbin;
     ATH_MSG_DEBUG("logEkin="<<logEkin<<" logEkin_nominal="<<logEkin_nominal<<" logEkin_next="<<logEkin_next<<" (rnd="<<rnd<<" < p(next)="<<numerator/denominator<<")? => orgbin="<<nextbin-1<<" selbin="<<bin);
   }
@@ -181,5 +183,3 @@ void TFCSParametrizationEkinSelectChain::unit_test(TFCSSimulationState* simulsta
   }  
   std::cout<<"==================================="<<std::endl<<std::endl;
 }
-
-
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSSimulationState.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSSimulationState.cxx
index fb7689e45f45c4e29bb29540e4f87290b272ecc0..36ba4ab873150795f5d65f57735144101e805493 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSSimulationState.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSSimulationState.cxx
@@ -2,6 +2,8 @@
   Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 */
 
+#include "CLHEP/Random/RandomEngine.h"
+
 #include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
 #include <iostream>
 
@@ -9,7 +11,8 @@
 //======= TFCSSimulationState =========
 //=============================================
 
-TFCSSimulationState::TFCSSimulationState()
+TFCSSimulationState::TFCSSimulationState(CLHEP::HepRandomEngine* randomEngine)
+  : m_randomEngine(randomEngine)
 {
   clear();
 }
@@ -47,5 +50,3 @@ void TFCSSimulationState::Print(Option_t *) const
     std::cout<<"  E"<<i<<"("<<CaloSampling::getSamplingName(i)<<")="<<E(i)<<" E"<<i<<"/E="<<Efrac(i)<<std::endl;
   }
 }
-
-
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSSimpleLateralShapeParametrization.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSSimpleLateralShapeParametrization.h
index 9b2ff0a30540fb3179dd3c62080e4561047b0fea..c1d83d13d6dc59499286c6631762da8ec622f6b8 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSSimpleLateralShapeParametrization.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSSimpleLateralShapeParametrization.h
@@ -10,12 +10,14 @@
 #include "TFile.h"
 #include "TH2F.h"
 #include "TF1.h"
-#include "TRandom3.h"
+
+namespace CLHEP {
+  class HepRandomEngine;
+}
 
 class TFCSSimpleLateralShapeParametrization:public TFCSLateralShapeParametrizationHitBase {
 public:
   TFCSSimpleLateralShapeParametrization(const char* name=nullptr, const char* title=nullptr);
-  ~TFCSSimpleLateralShapeParametrization();
 
   // simulated one hit position with weight that should be put into simulstate
   // sometime later all hit weights should be resacled such that their final sum is simulstate->E(sample)
@@ -27,7 +29,7 @@ public:
 
   bool Initialize(float input_sigma_x, float input_sigma_y);
 
-  void getHitXY(double &x, double &y);
+  void getHitXY(CLHEP::HepRandomEngine *engine, double &x, double &y);
 
   float getSigma_x(){return m_sigmaX;};
   float getSigma_y(){return m_sigmaY;};
@@ -37,14 +39,6 @@ private:
   float m_sigmaX;
   float m_sigmaY;
 
-  //float sigma2_x;
-  //float sigma2_y;
-
-  //float gaus_ratio;
-
-  TRandom3 *m_rnd;
-
-
   ClassDefOverride(TFCSSimpleLateralShapeParametrization,1)  //TFCSSimpleLateralShapeParametrization
 };
 
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSSimpleLateralShapeParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSSimpleLateralShapeParametrization.cxx
index fcd88efcdcc8ddc6b28548bd977dd4e923b6de3c..211febc6c967b5fd490388b09ba933dc57428402 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSSimpleLateralShapeParametrization.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSSimpleLateralShapeParametrization.cxx
@@ -2,8 +2,11 @@
   Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
 */
 
+#include "CLHEP/Random/RandGauss.h"
+
 #include "ISF_FastCaloSimParametrization/TFCSSimpleLateralShapeParametrization.h"
 #include "ISF_FastCaloSimEvent/TFCSExtrapolationState.h"
+#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
 #include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h"
 
 #include "TMath.h"
@@ -13,28 +16,26 @@
 //=============================================
 
 TFCSSimpleLateralShapeParametrization::TFCSSimpleLateralShapeParametrization(const char* name, const char* title) :
-  TFCSLateralShapeParametrizationHitBase(name,title),
-  m_rnd(0)
+  TFCSLateralShapeParametrizationHitBase(name,title)
 {
     m_sigmaX = 0;
     m_sigmaY = 0;
 }
 
-TFCSSimpleLateralShapeParametrization::~TFCSSimpleLateralShapeParametrization()
-{
-  if(m_rnd) delete m_rnd;
-}
-
 
-FCSReturnCode TFCSSimpleLateralShapeParametrization::simulate_hit(Hit& hit,TFCSSimulationState& /*simulstate*/,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* extrapol)
+FCSReturnCode TFCSSimpleLateralShapeParametrization::simulate_hit(Hit & hit, TFCSSimulationState &simulstate, const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* extrapol)
 {
+  if (!simulstate.randomEngine()) {
+    return FCSFatal;
+  }
+
   int cs=calosample();
   hit.eta()=0.5*( extrapol->eta(cs, CaloSubPos::SUBPOS_ENT) + extrapol->eta(cs, CaloSubPos::SUBPOS_EXT) );
   hit.phi()=0.5*( extrapol->phi(cs, CaloSubPos::SUBPOS_ENT) + extrapol->phi(cs, CaloSubPos::SUBPOS_EXT) );
   hit.E()*=1;
 
   double x, y;
-  getHitXY(x, y);
+  getHitXY(simulstate.randomEngine(), x, y);
 
   // delta_eta and delta_phi;
   double delta_eta = x;
@@ -49,10 +50,6 @@ FCSReturnCode TFCSSimpleLateralShapeParametrization::simulate_hit(Hit& hit,TFCSS
 
 bool TFCSSimpleLateralShapeParametrization::Initialize(float input_sigma_x, float input_sigma_y)
 {
-  // Setup random numbers
-  m_rnd = new TRandom3();
-  m_rnd->SetSeed(0);
-
   m_sigmaX = input_sigma_x;
   m_sigmaY = input_sigma_y;
   return true;
@@ -60,11 +57,6 @@ bool TFCSSimpleLateralShapeParametrization::Initialize(float input_sigma_x, floa
 
 bool TFCSSimpleLateralShapeParametrization::Initialize(const char* filepath, const char* histname)
 {
-    // Setup random numbers
-    m_rnd = new TRandom3();
-    m_rnd->SetSeed(0);
-
-
     // input file with histogram to fit
     TFile *f = new TFile(filepath);
     if (f == NULL) return false;
@@ -132,11 +124,8 @@ bool TFCSSimpleLateralShapeParametrization::Initialize(const char* filepath, con
     return true;
 }
 
-void TFCSSimpleLateralShapeParametrization::getHitXY(double &x, double &y)
+void TFCSSimpleLateralShapeParametrization::getHitXY(CLHEP::HepRandomEngine *engine, double &x, double &y)
 {
-
-
-    x = m_rnd->Gaus(0, m_sigmaX);
-    y = m_rnd->Gaus(0, m_sigmaY);
-
+    x = CLHEP::RandGauss::shoot(engine, 0, m_sigmaX);
+    y = CLHEP::RandGauss::shoot(engine, 0, m_sigmaY);
 }
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryLookup.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryLookup.cxx
index cb0bc117813c13fe8475b4c278d0aa2ff01d3ba3..c9d46ea4afcdf5148a8dedae4db1fc7153030da7 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryLookup.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/src/CaloGeometryLookup.cxx
@@ -6,7 +6,6 @@
 #include "ISF_FastCaloSimParametrization/CaloGeometryLookup.h"
 #include <TTree.h>
 #include <TVector2.h>
-#include <TRandom.h>
 #include <TCanvas.h>
 #include <TH2D.h>
 #include <TGraphErrors.h>
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/FastCaloSimSvcV2.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/FastCaloSimSvcV2.cxx
index b5c016701657065a24b3404923578f947e8bb660..23a02e310d80db83b9223e556cf1e24e97c2fa63 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/FastCaloSimSvcV2.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimServices/src/FastCaloSimSvcV2.cxx
@@ -28,10 +28,6 @@
 #include "StoreGate/StoreGateSvc.h"
 #include "StoreGate/StoreGate.h"
 
-
-#include "CLHEP/Random/RandFlat.h"
-#include "CLHEP/Random/RandGauss.h"
-
 #include "CaloEvent/CaloCellContainer.h"
 #include "CaloDetDescr/CaloDetDescrElement.h"
 #include "CaloDetDescr/CaloDetDescrManager.h"
@@ -206,7 +202,7 @@ StatusCode ISF::FastCaloSimSvcV2::simulate(const ISF::ISFParticle& isfp)
 
   TFCSExtrapolationState extrapol;
   m_FastCaloSimCaloExtrapolation->extrapolate(extrapol,&truth);
-  TFCSSimulationState simulstate;
+  TFCSSimulationState simulstate(m_randomEngine);
 
   FCSReturnCode status = m_param->simulate(simulstate, &truth, &extrapol);
   if (status != FCSSuccess) {
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PDFcreator.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PDFcreator.cxx
index 3f3b0ee0faa397aa2294e7277e4a13d975caa018..c8fb849fdc5ece6e9b6a89ada2eccf405986ced2 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PDFcreator.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PDFcreator.cxx
@@ -16,7 +16,6 @@
 #include "TFile.h"
 #include "TH2F.h"
 #include "TAxis.h"
-#include "TRandom.h"
 
 /*=========================================================================
  *  DESCRIPTION OF FUNCTION:
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PunchThroughTool.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PunchThroughTool.cxx
index e3c3d6464ed5efa3abc79039462f6beeaab5502d..a9cd1be6001a7b597a2e5891930327c5c6048df9 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PunchThroughTool.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PunchThroughTool.cxx
@@ -40,7 +40,6 @@
 #include "TH2F.h"
 #include "TAxis.h"
 #include "TH1.h"
-#include "TRandom.h"
 #include "TMath.h"
 
 // PathResolver