From 958818b70cf60afd94a94b0f84ea5fe15d33a1c2 Mon Sep 17 00:00:00 2001
From: Michael Duehrssen <michael.duehrssen@cern.ch>
Date: Tue, 16 Jan 2018 09:45:18 +0100
Subject: [PATCH] Update files for shape simulation and validation

Former-commit-id: a0f3efc34a5c0c3bbd5f62ff7247eb584c878cfa
---
 .../ISF_FastCaloSimEvent/TFCS1DFunction.h     |  4 +++
 .../TFCS1DFunctionHistogram.h                 |  4 +++
 .../TFCS1DFunctionRegression.h                |  4 +++
 .../TFCSEnergyParametrization.h               |  4 +++
 .../TFCSPCAEnergyParametrization.h            |  5 +++
 .../TFCSParametrization.h                     |  1 +
 .../TFCSParametrizationBase.h                 | 10 ++++--
 .../TFCSSimulationState.h                     | 29 +++++----------
 .../src/TFCSPCAEnergyParametrization.cxx      |  4 +++
 .../src/TFCSParametrization.cxx               |  5 +++
 .../src/TFCSSimulationState.cxx               | 20 ++++++++---
 .../CMakeLists.txt                            |  2 ++
 .../MeanAndRMS.h                              | 11 +++---
 .../TFCSHistoLateralShapeParametrization.h    |  7 ++--
 .../TFCSHitCellMappingWiggleEMB.h             |  2 +-
 .../TFCSLateralShapeParametrizationHitBase.h  |  2 +-
 .../TFCSLateralShapeParametrizationHitChain.h |  1 +
 .../Root/LinkDef.h                            |  2 ++
 .../TFCSHistoLateralShapeParametrization.cxx  | 25 +++++++++----
 .../Root/TFCSHitCellMapping.cxx               |  1 +
 .../Root/TFCSHitCellMappingWiggleEMB.cxx      |  2 ++
 .../Root/TFCSLateralShapeParametrization.cxx  | 12 +++----
 ...TFCSLateralShapeParametrizationHitBase.cxx |  6 +++-
 ...FCSLateralShapeParametrizationHitChain.cxx | 36 ++++++++++++++++---
 24 files changed, 145 insertions(+), 54 deletions(-)

diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunction.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunction.h
index a3b3e67ef4ca9..5ca5ceee51555 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunction.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunction.h
@@ -26,4 +26,8 @@ class TFCS1DFunction:public TObject
 
 };
 
+#if defined(__MAKECINT__)
+#pragma link C++ class TFCS1DFunction+;
+#endif
+
 #endif
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 b22154d0d0815..0bee476e5bedc 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunctionHistogram.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunctionHistogram.h
@@ -42,4 +42,8 @@ class TFCS1DFunctionHistogram:public TFCS1DFunction
 
 };
 
+#if defined(__MAKECINT__)
+#pragma link C++ class TFCS1DFunctionHistogram+;
+#endif
+
 #endif
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunctionRegression.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunctionRegression.h
index cbe3a9972579d..180fb26756164 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunctionRegression.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS1DFunctionRegression.h
@@ -34,4 +34,8 @@ class TFCS1DFunctionRegression:public TFCS1DFunction
 
 };
 
+#if defined(__MAKECINT__)
+#pragma link C++ class TFCS1DFunctionRegression+;
+#endif
+
 #endif
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyParametrization.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyParametrization.h
index c19a72bedea0b..52952a566bd5b 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyParametrization.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSEnergyParametrization.h
@@ -22,4 +22,8 @@ private:
   ClassDef(TFCSEnergyParametrization,1)  //TFCSEnergyParametrization
 };
 
+#if defined(__MAKECINT__)
+#pragma link C++ class TFCSEnergyParametrization+;
+#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 4836a047eabc8..ceb9798e1f221 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h
@@ -22,6 +22,7 @@ class TFCSPCAEnergyParametrization:public TFCSEnergyParametrization
   virtual void simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol);
   
   int n_pcabins()        { return m_numberpcabins; };
+  virtual int n_bins() {return m_numberpcabins;};
   IntArray* get_layers() { return m_RelevantLayers; };
   
   void P2X(TVectorD*, TVectorD* , TMatrixD* , int, double* , double* , int);
@@ -46,4 +47,8 @@ class TFCSPCAEnergyParametrization:public TFCSEnergyParametrization
  
 };
 
+#if defined(__MAKECINT__)
+#pragma link C++ class TFCSPCAEnergyParametrization+;
+#endif
+
 #endif
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrization.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrization.h
index 93d65130f9045..0cd4fcae80fc3 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrization.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrization.h
@@ -24,6 +24,7 @@ public:
   double eta_max() const {return m_eta_max;};
 
   void set_pdgid(int id);
+  void set_pdgid(const std::set< int > &ids);
   void add_pdgid(int id);
   void clear_pdgid();
 
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrizationBase.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrizationBase.h
index 2549ef81ae0cb..3a608757de82c 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrizationBase.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrizationBase.h
@@ -7,9 +7,11 @@
 
 #include <TNamed.h>
 #include <set>
-#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
-#include "ISF_FastCaloSimEvent/TFCSTruthState.h"
-#include "ISF_FastCaloSimEvent/TFCSExtrapolationState.h"
+
+class CaloGeometry;
+class TFCSSimulationState;
+class TFCSTruthState;
+class TFCSExtrapolationState;
 
 class TFCSParametrizationBase:public TNamed {
 public:
@@ -30,6 +32,8 @@ public:
   virtual double eta_min() const {return 100;};
   virtual double eta_max() const {return 100;};
 
+  virtual void set_geometry(CaloGeometry*) {};
+
   // Do some simulation. Result should be returned in simulstate
   // Simulate all energies in calo layers for energy parametrizations
   // Simulate one HIT for later shape parametrizations (TO BE DISCUSSED!)
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 d0e32ad8c8175..a341ddd88ceb5 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSSimulationState.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSSimulationState.h
@@ -7,8 +7,11 @@
 
 #include <TObject.h>
 #include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h"
+#include <map>
 #include <vector>
 
+class CaloDetDescrElement;
+
 class TFCSSimulationState:public TObject
 {
   public:
@@ -26,25 +29,10 @@ class TFCSSimulationState:public TObject
     void set_E(double E) { m_Etot=E; } ;
     void add_E(int sample,double Esample) { m_E[sample]+=Esample;m_Etot+=Esample; };
 
-    //empty function so far
-    //not sure if we should keep the function here or rather write a deposit_cell function or similar
-    class t_hit
-    {
-      public:
-      t_hit():m_eta(0),m_phi(0),m_weight(0) {};
-      t_hit(float eta, float phi, float weight):m_eta(eta),m_phi(phi),m_weight(weight) {};
-
-      float& eta() {return m_eta;};
-      float& phi() {return m_phi;};
-      float& weight() {return m_weight;};
-
-      private:
-      float m_eta,m_phi,m_weight;
-    };
-    virtual void deposit_HIT(int sample,const t_hit& hit) {m_hits[sample].push_back(hit);};
-    virtual void deposit_HIT(int sample,float hit_eta,float hit_phi,float hit_weight) {m_hits[sample].emplace_back(hit_eta,hit_phi,hit_weight);};
+    typedef std::map<const CaloDetDescrElement*,float> t_cellmap;
 
-    std::vector< t_hit >& get_hits(int sample) {return m_hits[sample];};
+    t_cellmap& cells() {return m_cells;};
+    void deposit(const CaloDetDescrElement* cellele, float E);
     
     void Print(Option_t *option="") const;
 
@@ -52,11 +40,12 @@ class TFCSSimulationState:public TObject
   private:
     int    m_Ebin;
     double m_Etot;
+    // TO BE CLEANED UP! SHOULD ONLY STORE EITHER E OR EFRAC!!!
     double m_E[CaloCell_ID_FCS::MaxSample];
     double m_Efrac[CaloCell_ID_FCS::MaxSample];
     
-    std::vector< std::vector< t_hit > > m_hits;
-
+    t_cellmap m_cells;
+    
   ClassDef(TFCSSimulationState,1)  //TFCSSimulationState
 };
 
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
index 20098f421d51a..5ef293a98156d 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSPCAEnergyParametrization.cxx
@@ -5,6 +5,9 @@
 #include "ISF_FastCaloSimEvent/TFCSPCAEnergyParametrization.h"
 #include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h"
 
+#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
+#include "ISF_FastCaloSimEvent/TFCSExtrapolationState.h"
+
 #include "TFile.h"
 #include <iostream>
 #include "TKey.h"
@@ -12,6 +15,7 @@
 #include "TRandom3.h"
 #include "TMatrixD.h"
 #include "TMatrixDSymEigen.h"
+#include "TMath.h"
 
 //=============================================
 //======= TFCSPCAEnergyParametrization =========
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrization.cxx
index 6d4caf53fe878..d4e3bb1bfc324 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrization.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSParametrization.cxx
@@ -18,6 +18,11 @@ void TFCSParametrization::set_pdgid(int id)
   m_pdgid.insert(id);
 }
 
+void TFCSParametrization::set_pdgid(const std::set< int > &ids)
+{
+  m_pdgid=ids;
+}
+
 void TFCSParametrization::add_pdgid(int id)
 {
   m_pdgid.insert(id);
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSSimulationState.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSSimulationState.cxx
index 6327f9150dc13..fb7689e45f45c 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSSimulationState.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSSimulationState.cxx
@@ -18,21 +18,33 @@ void TFCSSimulationState::clear()
 {
   m_Ebin=-1;
   m_Etot=0;
-  m_hits.resize(CaloCell_ID_FCS::MaxSample);
   for(int i=0;i<CaloCell_ID_FCS::MaxSample;++i)
   {
     m_E[i]=0;
     m_Efrac[i]=0;
-    m_hits[i].clear();
   }
 }
 
+void TFCSSimulationState::deposit(const CaloDetDescrElement* cellele, float E)
+{
+  //std::cout<<"TFCSSimulationState::deposit: cellele="<<cellele<<" E="<<E;
+  auto mele=m_cells.find(cellele);
+  if(mele==m_cells.end()) {
+    m_cells.emplace(cellele,0);
+    mele=m_cells.find(cellele);
+  }  
+  //std::cout<<" Ebefore="<<mele->second;
+  m_cells[cellele]+=E;
+  //std::cout<<" Eafter="<<mele->second;
+  //std::cout<<std::endl;
+}
+
 void TFCSSimulationState::Print(Option_t *) const
 {
-  std::cout<<"Ebin="<<m_Ebin<<" E="<<E()<<std::endl;
+  std::cout<<"Ebin="<<m_Ebin<<" E="<<E()<<" #cells="<<m_cells.size()<<std::endl;
   for(int i=0;i<CaloCell_ID_FCS::MaxSample;++i) if(E(i)!=0)
   {
-    std::cout<<"  E"<<i<<"("<<CaloSampling::getSamplingName(i)<<")="<<E(i)<<" E"<<i<<"/E="<<Efrac(i)<<" #hits="<<m_hits[i].size()<<std::endl;
+    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/CMakeLists.txt b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/CMakeLists.txt
index feaaef3dbe3d1..a46b78db6f02c 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/CMakeLists.txt
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/CMakeLists.txt
@@ -60,6 +60,8 @@ atlas_add_root_dictionary( ISF_FastCaloSimParametrizationLib
                            ISF_FastCaloSimParametrization/TFCSNNLateralShapeParametrization.h
                            ISF_FastCaloSimParametrization/TFCSSimpleLateralShapeParametrization.h
                            ISF_FastCaloSimParametrization/TFCSHistoLateralShapeParametrization.h
+                           ISF_FastCaloSimParametrization/TFCSLateralShapeParametrizationHitChain.h
+                           ISF_FastCaloSimParametrization/TFCSLateralShapeParametrizationHitBase.h
                            ISF_FastCaloSimParametrization/TreeReader.h
                            ISF_FastCaloSimParametrization/FCS_Cell.h
                            Root/LinkDef.h
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/MeanAndRMS.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/MeanAndRMS.h
index f9e28f262d66e..854c774e3cdd1 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/MeanAndRMS.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/MeanAndRMS.h
@@ -10,18 +10,21 @@
 class MeanAndRMS {
 public :
   MeanAndRMS():m_w(0),m_wx(0),m_wx2(0) {};
-  virtual ~MeanAndRMS() {};
+  MeanAndRMS(const double xadd, const double weight=1):m_w(weight),m_wx(weight*xadd),m_wx2(weight*xadd*xadd) {};
+  MeanAndRMS(const MeanAndRMS& ref):m_w(ref.m_w),m_wx(ref.m_wx),m_wx2(ref.m_wx2) {};
    
   MeanAndRMS& add(double xadd,double weight=1) {m_wx+=weight*xadd;m_wx2+=weight*xadd*xadd;m_w+=weight;return *this;};
   MeanAndRMS& operator+=(double xadd) {return add(xadd);};
   MeanAndRMS& operator-=(double xadd) {return add(-xadd);};
+  MeanAndRMS& operator=(const MeanAndRMS& ref) {m_w=ref.m_w;m_wx=ref.m_wx;m_wx2=ref.m_wx2;return *this;};
   
   double sum_weight() const {return m_w;};
   double mean()       const {if(m_w!=0) return m_wx/m_w; else return 0;};
-  double rms2()       const {double x=mean();return m_wx2/m_w - x*x;};
+  double mean2()      const {double x=mean();return x*x;};
+  double rms2()       const {if(m_w!=0) return m_wx2/m_w - mean2(); else return 0;};
   double rms()        const {double r2=rms2();if(r2>=0) return sqrt(r2); else return 0;};
-  double mean_error() const {return rms()/sqrt(m_w);};
-  double rms_error()  const {return rms()/sqrt(2*m_w);};
+  double mean_error() const {if(m_w>0) return rms()/sqrt(m_w); else return 0;};
+  double rms_error()  const {if(m_w>0) return rms()/sqrt(2*m_w); else return 0;};
 
   operator double() const { return mean(); }
 protected:
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSHistoLateralShapeParametrization.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSHistoLateralShapeParametrization.h
index d8256e20fd033..694d090870ee2 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSHistoLateralShapeParametrization.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSHistoLateralShapeParametrization.h
@@ -5,20 +5,21 @@
 #ifndef TFCSHistoLateralShapeParametrization_h
 #define TFCSHistoLateralShapeParametrization_h
 
-#include "ISF_FastCaloSimParametrization/TFCSLateralShapeParametrization.h"
+#include "ISF_FastCaloSimParametrization/TFCSLateralShapeParametrizationHitBase.h"
 
 #include "TH2.h"
 #include "TRandom3.h"
 
 
-class TFCSHistoLateralShapeParametrization:public TFCSLateralShapeParametrization {
+class TFCSHistoLateralShapeParametrization:public TFCSLateralShapeParametrizationHitBase {
 public:
   TFCSHistoLateralShapeParametrization(const char* name=0, const char* title=0);
 
+  int get_number_of_hits(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const;
   // 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)
   // someone also needs to map all hits into cells
-  virtual void simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol);
+  virtual void simulate_hit(t_hit& hit,TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol);
 
   // Init and fill sigma
   bool Initialize(TH2* hist);
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSHitCellMappingWiggleEMB.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSHitCellMappingWiggleEMB.h
index f624e13359743..f3c7a2f0e5ffb 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSHitCellMappingWiggleEMB.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSHitCellMappingWiggleEMB.h
@@ -13,7 +13,7 @@ class TFCSHitCellMappingWiggleEMB:public TFCSLateralShapeParametrizationHitBase
 public:
   TFCSHitCellMappingWiggleEMB(const char* name=0, const char* title=0, CaloGeometry* geo=0);
   
-  void set_geometry(CaloGeometry* geo) {m_geo=geo;};
+  virtual void set_geometry(CaloGeometry* geo) {m_geo=geo;};
   CaloGeometry* get_geometry() {return m_geo;};
 
   // simulated one hit position with weight that should be put into simulstate
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSLateralShapeParametrizationHitBase.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSLateralShapeParametrizationHitBase.h
index 00c05ddf183d4..14afccbb2d651 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSLateralShapeParametrizationHitBase.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSLateralShapeParametrizationHitBase.h
@@ -30,7 +30,7 @@ public:
   class t_hit
   {
     public:
-    t_hit():m_eta(0),m_phi(0),m_E(0) {};
+    t_hit():m_eta(0),m_phi(0),m_E(0) {}; // for hits with the same energy, m_E should normalized to E(layer)/nhit
     t_hit(float eta, float phi, float E):m_eta(eta),m_phi(phi),m_E(E) {};
 
     float& eta() {return m_eta;};
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSLateralShapeParametrizationHitChain.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSLateralShapeParametrizationHitChain.h
index 2c0c68ada2fb5..656c2d1f9f751 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSLateralShapeParametrizationHitChain.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/ISF_FastCaloSimParametrization/TFCSLateralShapeParametrizationHitChain.h
@@ -13,6 +13,7 @@ class TFCSLateralShapeParametrizationHitBase;
 class TFCSLateralShapeParametrizationHitChain:public TFCSLateralShapeParametrization {
 public:
   TFCSLateralShapeParametrizationHitChain(const char* name=0, const char* title=0);
+  TFCSLateralShapeParametrizationHitChain(TFCSLateralShapeParametrizationHitBase* hitsim);
 
   virtual void simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol);
 
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/LinkDef.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/LinkDef.h
index 0874ea551a23f..22a02326693eb 100755
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/LinkDef.h
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/LinkDef.h
@@ -7,6 +7,8 @@
 #pragma link C++ class TFCSNNLateralShapeParametrization+;
 #pragma link C++ class TFCSSimpleLateralShapeParametrization+;
 #pragma link C++ class TFCSHistoLateralShapeParametrization+;
+#pragma link C++ class TFCSLateralShapeParametrizationHitChain.h+;
+#pragma link C++ class TFCSLateralShapeParametrizationHitBase.h+;
 #pragma link C++ class MeanAndRMS;
 #pragma link C++ class TreeReader;
 #pragma link C++ class EnergyParametrizationValidation;
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSHistoLateralShapeParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSHistoLateralShapeParametrization.cxx
index 2835b51db388e..3d9d365008524 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSHistoLateralShapeParametrization.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSHistoLateralShapeParametrization.cxx
@@ -4,31 +4,41 @@
 
 #include "ISF_FastCaloSimParametrization/TFCSHistoLateralShapeParametrization.h"
 #include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h"
+#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
+#include "ISF_FastCaloSimEvent/TFCSExtrapolationState.h"
 
 #include "TFile.h"
+#include "TMath.h"
+
+#include <iostream>
 
 //=============================================
 //======= TFCSHistoLateralShapeParametrization =========
 //=============================================
 
 TFCSHistoLateralShapeParametrization::TFCSHistoLateralShapeParametrization(const char* name, const char* title) :
-  TFCSLateralShapeParametrization(name,title),
+  TFCSLateralShapeParametrizationHitBase(name,title),
   m_hist(0),m_rnd(0)
 {
 }
 
-void TFCSHistoLateralShapeParametrization::simulate(TFCSSimulationState& simulstate,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* extrapol)
+int TFCSHistoLateralShapeParametrization::get_number_of_hits(TFCSSimulationState& simulstate,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* /*extrapol*/) const
+{
+  return gRandom->Poisson(m_hist->Integral());
+}
+
+void TFCSHistoLateralShapeParametrization::simulate_hit(t_hit& hit,TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol)
 {
   int cs=calosample();
   double center_eta=0.5*( extrapol->eta(cs, CaloSubPos::SUBPOS_ENT) + extrapol->eta(cs, CaloSubPos::SUBPOS_EXT) );
   double center_phi=0.5*( extrapol->phi(cs, CaloSubPos::SUBPOS_ENT) + extrapol->phi(cs, CaloSubPos::SUBPOS_EXT) );
   double center_r=0.5*( extrapol->r(cs, CaloSubPos::SUBPOS_ENT) + extrapol->r(cs, CaloSubPos::SUBPOS_EXT) );
   double center_z=0.5*( extrapol->z(cs, CaloSubPos::SUBPOS_ENT) + extrapol->z(cs, CaloSubPos::SUBPOS_EXT) );
-  double hit_weight=1;
+  hit.E()*=1;
 
   double alpha, r;
   
-  m_hist->GetRandom2(alpha, r);
+  m_hist->GetRandom2(r,alpha);
   double delta_eta_mm = r * cos(alpha);
   double delta_phi_mm = r * sin(alpha);
 
@@ -39,9 +49,10 @@ void TFCSHistoLateralShapeParametrization::simulate(TFCSSimulationState& simulst
   double delta_eta = delta_eta_mm / eta_jakobi / dist000;
   double delta_phi = delta_phi_mm / center_r;
 
-  double hit_eta = center_eta + delta_eta;
-  double hit_phi = center_phi + delta_phi;
-  simulstate.deposit_HIT(cs,hit_eta,hit_phi,hit_weight);
+  hit.eta() = center_eta + delta_eta;
+  hit.phi() = center_phi + delta_phi;
+  //simulstate.deposit_HIT(cs,hit_eta,hit_phi,hit_weight);
+  //std::cout<<"TFCSHistoLateralShapeParametrization::simulate_hit: E="<<hit.E()<<" cs="<<cs<<" eta="<<hit.eta()<<" phi="<<hit.phi()<<" r="<<r<<" alpha="<<alpha<<" delta_eta_mm="<<delta_eta_mm<<" delta_phi_mm="<<delta_phi_mm<<std::endl;
 }
 
 
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSHitCellMapping.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSHitCellMapping.cxx
index e4591698ae344..a5ebf3d47424c 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSHitCellMapping.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSHitCellMapping.cxx
@@ -4,6 +4,7 @@
 
 #include "ISF_FastCaloSimParametrization/TFCSHitCellMapping.h"
 #include "ISF_FastCaloSimParametrization/CaloGeometry.h"
+#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
 #include <iostream>
 
 //=============================================
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSHitCellMappingWiggleEMB.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSHitCellMappingWiggleEMB.cxx
index 9b232a56487ff..d48a9d2ee0df4 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSHitCellMappingWiggleEMB.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSHitCellMappingWiggleEMB.cxx
@@ -4,8 +4,10 @@
 
 #include "ISF_FastCaloSimParametrization/TFCSHitCellMappingWiggleEMB.h"
 #include "ISF_FastCaloSimParametrization/CaloGeometry.h"
+#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
 #include <iostream>
 #include "TRandom3.h"
+#include "TVector2.h"
 
 //=============================================
 //======= TFCSHitCellMappingWiggleEMB =========
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSLateralShapeParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSLateralShapeParametrization.cxx
index 783fd7c04f112..037bb1979c718 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSLateralShapeParametrization.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSLateralShapeParametrization.cxx
@@ -24,14 +24,14 @@ void TFCSLateralShapeParametrization::set_calosample(int cs)
   m_calosample=cs;
 }
 
-void TFCSLateralShapeParametrization::simulate(TFCSSimulationState& simulstate,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* extrapol)
+void TFCSLateralShapeParametrization::simulate(TFCSSimulationState& /*simulstate*/,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* /*extrapol*/)
 {
-  int cs=calosample();
-  double hit_eta=0.5*( extrapol->eta(cs, CaloSubPos::SUBPOS_ENT) + extrapol->eta(cs, CaloSubPos::SUBPOS_EXT) );
-  double hit_phi=0.5*( extrapol->phi(cs, CaloSubPos::SUBPOS_ENT) + extrapol->phi(cs, CaloSubPos::SUBPOS_EXT) );
-  double hit_weight=1;
+//  int cs=calosample();
+//  double hit_eta=0.5*( extrapol->eta(cs, CaloSubPos::SUBPOS_ENT) + extrapol->eta(cs, CaloSubPos::SUBPOS_EXT) );
+//  double hit_phi=0.5*( extrapol->phi(cs, CaloSubPos::SUBPOS_ENT) + extrapol->phi(cs, CaloSubPos::SUBPOS_EXT) );
+//  double hit_weight=1;
 
-  simulstate.deposit_HIT(cs,hit_eta,hit_phi,hit_weight);
+//  simulstate.deposit_HIT(cs,hit_eta,hit_phi,hit_weight);
 }
 
 void TFCSLateralShapeParametrization::Print(Option_t *option) const
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSLateralShapeParametrizationHitBase.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSLateralShapeParametrizationHitBase.cxx
index 07fa75b478b75..6cb4ea49b47aa 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSLateralShapeParametrizationHitBase.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSLateralShapeParametrizationHitBase.cxx
@@ -4,6 +4,10 @@
 
 #include "ISF_FastCaloSimParametrization/TFCSLateralShapeParametrizationHitBase.h"
 #include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h"
+
+#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
+#include "ISF_FastCaloSimEvent/TFCSExtrapolationState.h"
+
 #include <iostream>
 
 //=============================================
@@ -34,7 +38,7 @@ void TFCSLateralShapeParametrizationHitBase::simulate_hit(t_hit& hit,TFCSSimulat
   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;
+  hit.E()*=1;
 }
 
 void TFCSLateralShapeParametrizationHitBase::Print(Option_t *option) const
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSLateralShapeParametrizationHitChain.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSLateralShapeParametrizationHitChain.cxx
index b311297f0ab2d..ab4991234f3f1 100644
--- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSLateralShapeParametrizationHitChain.cxx
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/Root/TFCSLateralShapeParametrizationHitChain.cxx
@@ -7,6 +7,8 @@
 #include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h"
 #include <iostream>
 
+#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
+
 //=============================================
 //======= TFCSLateralShapeParametrization =========
 //=============================================
@@ -15,11 +17,31 @@ TFCSLateralShapeParametrizationHitChain::TFCSLateralShapeParametrizationHitChain
 {
 }
 
+TFCSLateralShapeParametrizationHitChain::TFCSLateralShapeParametrizationHitChain(TFCSLateralShapeParametrizationHitBase* hitsim):TFCSLateralShapeParametrization(TString("hit_chain_")+hitsim->GetName(),TString("hit chain for ")+hitsim->GetTitle())
+{
+  set_calosample(hitsim->calosample());
+  
+  set_Ekin_bin(hitsim->Ekin_bin());
+  set_Ekin_nominal(hitsim->Ekin_nominal());
+  set_Ekin_min(hitsim->Ekin_min());
+  set_Ekin_max(hitsim->Ekin_max());
+  
+  set_eta_nominal(hitsim->eta_nominal());
+  set_eta_min(hitsim->eta_min());
+  set_eta_max(hitsim->eta_max());
+
+  set_pdgid(hitsim->pdgid());
+
+  m_chain.push_back(hitsim);
+}
+
 void TFCSLateralShapeParametrizationHitChain::set_geometry(CaloGeometry* geo)
 {
   for(TFCSLateralShapeParametrizationHitBase* hitsim : m_chain) hitsim->set_geometry(geo);
 }
 
+
+
 int TFCSLateralShapeParametrizationHitChain::get_number_of_hits(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const
 {
   for(TFCSLateralShapeParametrizationHitBase* hitsim : m_chain) {
@@ -32,10 +54,16 @@ int TFCSLateralShapeParametrizationHitChain::get_number_of_hits(TFCSSimulationSt
 
 void TFCSLateralShapeParametrizationHitChain::simulate(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol)
 {
-  TFCSLateralShapeParametrizationHitBase::t_hit hit;
-  for(TFCSLateralShapeParametrizationHitBase* hitsim : m_chain) {
-    hitsim->simulate_hit(hit,simulstate,truth,extrapol);
-  } 
+  // Call get_number_of_hits() only once, as it could contain a random number
+  int nhit=get_number_of_hits(simulstate,truth,extrapol);
+  float Ehit=simulstate.E(calosample())/nhit;
+  for(int i=0;i<nhit;++i) {
+    TFCSLateralShapeParametrizationHitBase::t_hit hit; 
+    hit.E()=Ehit;
+    for(TFCSLateralShapeParametrizationHitBase* hitsim : m_chain) {
+      hitsim->simulate_hit(hit,simulstate,truth,extrapol);
+    } 
+  }  
 }
 
 void TFCSLateralShapeParametrizationHitChain::Print(Option_t *option) const
-- 
GitLab