Add functionality for weighted hits in the lateral shape simulation

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt

@@ -63,6 +63,8 @@ atlas_add_root_dictionary( ISF_FastCaloSimEvent _dictSource
                            ISF_FastCaloSimEvent/TFCSCenterPositionCalculation.h
                            ISF_FastCaloSimEvent/TFCSHistoLateralShapeParametrization.h
                            ISF_FastCaloSimEvent/TFCSHistoLateralShapeParametrizationFCal.h
+                           ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h
+                           ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeightHitAndMiss.h
                            ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitNumberFromE.h
                            ISF_FastCaloSimEvent/TFCSHitCellMapping.h
                            ISF_FastCaloSimEvent/TFCSHitCellMappingFCal.h

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/LinkDef.h

@@ -42,6 +42,8 @@
 #include "ISF_FastCaloSimEvent/TFCSCenterPositionCalculation.h"
 #include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeParametrization.h"
 #include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeParametrizationFCal.h"
+#include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h"
+#include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeightHitAndMiss.h"
 #include "ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitNumberFromE.h"
 #include "ISF_FastCaloSimEvent/TFCSHitCellMapping.h"
 #include "ISF_FastCaloSimEvent/TFCSHitCellMappingFCal.h"
@@ -150,6 +152,8 @@
 #pragma link C++ class TFCSCenterPositionCalculation+;
 #pragma link C++ class TFCSHistoLateralShapeParametrization+;
 #pragma link C++ class TFCSHistoLateralShapeParametrizationFCal+;
+#pragma link C++ class TFCSHistoLateralShapeWeight+;
+#pragma link C++ class TFCSHistoLateralShapeWeightHitAndMiss+;
 #pragma link C++ class TFCSLateralShapeParametrizationHitNumberFromE+;
 #pragma link C++ class TFCSHitCellMapping+;
 #pragma link C++ class TFCSHitCellMappingFCal+;

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h

+/*
+  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef TFCSHistoLateralShapeWeight_h
+#define TFCSHistoLateralShapeWeight_h
+
+#include "ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitBase.h"
+
+class TH1;
+
+class TFCSHistoLateralShapeWeight:public TFCSLateralShapeParametrizationHitBase {
+public:
+  TFCSHistoLateralShapeWeight(const char* name=nullptr, const char* title=nullptr);
+  ~TFCSHistoLateralShapeWeight();
+
+  /// weight the energy of one hit in order to generate fluctuations
+  virtual FCSReturnCode simulate_hit(Hit& hit,TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol);
+
+  /// Init from histogram. The integral of the histogram is used as number of expected hits to be generated
+  bool Initialize(TH1* hist);
+  
+  void Print(Option_t *option = "") const;
+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};
+
+  ClassDef(TFCSHistoLateralShapeWeight,1)  //TFCSHistoLateralShapeWeight
+};
+
+#if defined(__ROOTCLING__) && defined(__FastCaloSimStandAlone__)
+#pragma link C++ class TFCSHistoLateralShapeWeight+;
+#endif
+
+#endif

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeightHitAndMiss.h

+/*
+  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef TFCSHistoLateralShapeWeightHitAndMiss_h
+#define TFCSHistoLateralShapeWeightHitAndMiss_h
+
+#include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h"
+
+class TH1;
+
+class TFCSHistoLateralShapeWeightHitAndMiss:public TFCSHistoLateralShapeWeight {
+public:
+  TFCSHistoLateralShapeWeightHitAndMiss(const char* name=nullptr, const char* title=nullptr):TFCSHistoLateralShapeWeight(name,title) {};
+  ~TFCSHistoLateralShapeWeightHitAndMiss() {};
+
+  /// weight the energy of one hit in order to generate fluctuations. If the hit energy is 0, discard the hit
+  virtual FCSReturnCode simulate_hit(Hit& hit,TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol);
+
+private:
+
+  ClassDef(TFCSHistoLateralShapeWeightHitAndMiss,1)  //TFCSHistoLateralShapeWeightHitAndMiss
+};
+
+#if defined(__ROOTCLING__) && defined(__FastCaloSimStandAlone__)
+#pragma link C++ class TFCSHistoLateralShapeWeightHitAndMiss+;
+#endif
+
+#endif

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSHitCellMappingWiggle.h

@@ -26,7 +26,7 @@ public:
   inline double get_bin_low_edge(int bin) const {return m_bin_low_edge[bin];};
   inline double get_bin_up_edge(int bin) const {return m_bin_low_edge[bin+1];};
   
-  inline const TFCS1DFunction* get_function(int bin) {return m_functions[bin];};
+  inline const TFCS1DFunction* get_function(int bin) const {return m_functions[bin];};
   const std::vector< const TFCS1DFunction* > get_functions() {return m_functions;};
   const std::vector< float > get_bin_low_edges() {return m_bin_low_edge;};
 

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitBase.h

@@ -44,6 +44,10 @@ public:
       m_z=0.;
       m_E=0.;
       m_useXYZ=false;
+      m_center_r=0;
+      m_center_z=0;
+      m_center_eta=0;
+      m_center_phi=0;
     }
 
     inline float& eta() {return m_eta_x;};

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSParametrizationBase.h

@@ -141,12 +141,12 @@ public:
   static void DoCleanup();
 
 protected:
-  const double init_Ekin_nominal=0;
-  const double init_Ekin_min=0;
-  const double init_Ekin_max=14000000;
-  const double init_eta_nominal=0;
-  const double init_eta_min=-100;
-  const double init_eta_max=100;
+  static constexpr double init_Ekin_nominal=0;
+  static constexpr double init_Ekin_min=0;
+  static constexpr double init_Ekin_max=14000000;
+  static constexpr double init_eta_nominal=0;
+  static constexpr double init_eta_min=-100;
+  static constexpr double init_eta_max=100;
 
   static std::vector< TFCSParametrizationBase* > s_cleanup_list;
 
@@ -196,7 +196,7 @@ private:
 private:
   static std::set< int > s_no_pdgid;
 
-  ClassDef(TFCSParametrizationBase,1)  //TFCSParametrizationBase
+  ClassDef(TFCSParametrizationBase,2)  //TFCSParametrizationBase
 };
 
 #if defined(__ROOTCLING__) && defined(__FastCaloSimStandAlone__)

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCS1DFunction.cxx

@@ -113,7 +113,7 @@ void TFCS1DFunction::unit_test(TH1* hist,TFCS1DFunction* rtof,int nrnd,TH1* hist
     std::cout<<"rnd0="<<rnd[0]<<" -> x="<<value[0]<<std::endl;
   }
 
-  TH1* hist_val;
+  TH1* hist_val=nullptr;
   if(histfine) hist_val=(TH1*)histfine->Clone(TString(hist->GetName())+"hist_val");
    else hist_val=(TH1*)hist->Clone(TString(hist->GetName())+"hist_val");
   double weightfine=hist_val->Integral()/nrnd;
@@ -142,7 +142,6 @@ void TFCS1DFunction::unit_test(TH1* hist,TFCS1DFunction* rtof,int nrnd,TH1* hist
     float val=hist_diff->GetBinContent(ix);
     float err=hist_diff->GetBinError(ix);
     if(err>0) hist_pull->Fill(val/err);
-    //std::cout<<"x="<<hist->GetBinCenter(ix)<<" : pull val="<<val<<" err="<<err<<std::endl;
   }
   
 //Screen output in athena won't make sense and would require linking of additional libraries

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCS1DFunctionInt32Histogram.cxx

@@ -53,8 +53,7 @@ void TFCS1DFunctionInt32Histogram::Initialize(const TH1* hist)
   
   for(ibin=0;ibin<nbins;++ibin) {
     m_HistoContents[ibin]=s_MaxValue*(temp_HistoContents[ibin]/integral);
-    //std::cout<<"bin="<<ibin<<" val="<<m_HistoContents[ibin]<<std::endl;
-  }  
+  }
 }
 
 double TFCS1DFunctionInt32Histogram::rnd_to_fct(double rnd) const

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeWeight.cxx

+/*
+  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "CLHEP/Random/RandFlat.h"
+
+#include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h"
+#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
+
+#include "TH1.h"
+#include "TVector2.h"
+#include "TMath.h"
+
+//=============================================
+//======= TFCSHistoLateralShapeWeight =========
+//=============================================
+
+TFCSHistoLateralShapeWeight::TFCSHistoLateralShapeWeight(const char* name, const char* title) :
+  TFCSLateralShapeParametrizationHitBase(name,title)
+{
+}
+
+TFCSHistoLateralShapeWeight::~TFCSHistoLateralShapeWeight()
+{
+  if(m_hist) delete m_hist;
+}
+
+FCSReturnCode TFCSHistoLateralShapeWeight::simulate_hit(Hit& hit,TFCSSimulationState& simulstate,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* /*extrapol*/)
+{
+  if (!simulstate.randomEngine()) {
+    return FCSFatal;
+  }
+
+  const double center_eta = hit.center_eta(); 
+  const double center_phi = hit.center_phi();
+  const double center_r   = hit.center_r();
+  const double center_z   = hit.center_z();
+
+  const float dist000    = TMath::Sqrt(center_r * center_r + center_z * center_z);
+  const float eta_jakobi = TMath::Abs(2.0 * TMath::Exp(-center_eta) / (1.0 + TMath::Exp(-2 * center_eta)));
+
+  const float delta_eta = hit.eta()-center_eta;
+  const float delta_phi = hit.phi()-center_phi;
+  const float delta_eta_mm = delta_eta * eta_jakobi * dist000;
+  const float delta_phi_mm = delta_phi * center_r;
+  const float delta_r_mm = TMath::Sqrt(delta_eta_mm*delta_eta_mm+delta_phi_mm*delta_phi_mm);
+  
+  //TODO: delta_r_mm should perhaps be cached in hit
+
+  Int_t bin=m_hist->FindBin(delta_r_mm);
+  if(bin<1) bin=1;
+  if(bin>m_hist->GetNbinsX()) bin=m_hist->GetNbinsX();
+  float weight=m_hist->GetBinContent(bin);
+  hit.E()*=weight;
+
+  ATH_MSG_DEBUG("HIT: E="<<hit.E()<<" dR_mm="<<delta_r_mm<<" weight="<<weight);
+  return FCSSuccess;
+}
+
+
+bool TFCSHistoLateralShapeWeight::Initialize(TH1* hist)
+{
+  if(!hist) return false;
+  if(m_hist) delete m_hist;
+	m_hist=(TH1*)hist->Clone(TString("TFCSHistoLateralShapeWeight_")+hist->GetName());
+	m_hist->SetDirectory(0);
+
+  return true;
+}
+
+void TFCSHistoLateralShapeWeight::Print(Option_t *option) const
+{
+  TString opt(option);
+  bool shortprint=opt.Index("short")>=0;
+  bool longprint=msgLvl(MSG::DEBUG) || (msgLvl(MSG::INFO) && !shortprint);
+  TString optprint=opt;optprint.ReplaceAll("short","");
+  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()<<"]");
+     else ATH_MSG_INFO(optprint <<"  no Histogram");
+  }  
+}
+

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHistoLateralShapeWeightHitAndMiss.cxx

+/*
+  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "CLHEP/Random/RandFlat.h"
+
+#include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeightHitAndMiss.h"
+#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
+
+#include "TH1.h"
+#include "TVector2.h"
+#include "TMath.h"
+
+//=============================================
+//======= TFCSHistoLateralShapeWeightHitAndMiss =========
+//=============================================
+
+FCSReturnCode TFCSHistoLateralShapeWeightHitAndMiss::simulate_hit(Hit& hit,TFCSSimulationState& simulstate,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* /*extrapol*/)
+{
+  if (!simulstate.randomEngine()) {
+    return FCSFatal;
+  }
+
+  const double center_eta = hit.center_eta(); 
+  const double center_phi = hit.center_phi();
+  const double center_r   = hit.center_r();
+  const double center_z   = hit.center_z();
+
+  const float dist000    = TMath::Sqrt(center_r * center_r + center_z * center_z);
+  const float eta_jakobi = TMath::Abs(2.0 * TMath::Exp(-center_eta) / (1.0 + TMath::Exp(-2 * center_eta)));
+
+  const float delta_eta = hit.eta()-center_eta;
+  const float delta_phi = hit.phi()-center_phi;
+  const float delta_eta_mm = delta_eta * eta_jakobi * dist000;
+  const float delta_phi_mm = delta_phi * center_r;
+  const float delta_r_mm = TMath::Sqrt(delta_eta_mm*delta_eta_mm+delta_phi_mm*delta_phi_mm);
+  
+  //TODO: delta_r_mm should perhaps be cached in hit
+
+  Int_t bin=m_hist->FindBin(delta_r_mm);
+  if(bin<1) bin=1;
+  if(bin>m_hist->GetNbinsX()) bin=m_hist->GetNbinsX();
+  float weight=m_hist->GetBinContent(bin);
+  if(weight<=1) {
+    //if weight<=1, give lower energy to hit. 
+    //TFCSLateralShapeParametrizationHitChain needs to be able to generate more hits in this case
+    hit.E()*=weight;
+  } else {
+    //if weight>1, accept only 1/weight events, but give them a higher energy increased by weight. 
+    //This leads to larger fluctuations, while keeping the shape unchanged.
+    float prob=1.0/weight;
+    float rnd=CLHEP::RandFlat::shoot(simulstate.randomEngine());
+    if(rnd<prob) hit.E()*=weight;
+     else hit.E()=0;
+  }  
+
+  ATH_MSG_DEBUG("HIT: E="<<hit.E()<<" dR_mm="<<delta_r_mm<<" weight="<<weight);
+  return FCSSuccess;
+}
+
+

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSHitCellMappingWiggle.cxx

@@ -47,7 +47,7 @@ void TFCSHitCellMappingWiggle::initialize(TFCS1DFunction* func)
 void TFCSHitCellMappingWiggle::initialize(const std::vector< const TFCS1DFunction* >& functions, const std::vector< float >& bin_low_edges)
 {
   if(functions.size()+1!=bin_low_edges.size()) {
-    ATH_MSG_ERROR("Using "<<functions.size()<<" functions needs "<<functions.size()+1<<" bins, but got "<<bin_low_edges.size()<<"bins");
+    ATH_MSG_ERROR("Using "<<functions.size()<<" functions needs "<<functions.size()+1<<" bin low edges, but got "<<bin_low_edges.size()<<"bins");
     return;
   }
   for(auto function : m_functions) if(function) delete function;
@@ -125,9 +125,9 @@ void TFCSHitCellMappingWiggle::Print(Option_t *option) const
   TString optprint=opt;optprint.ReplaceAll("short","");
   
   if(longprint) {
-    ATH_MSG(INFO) << optprint <<"  "<<get_number_of_bins()<<" functions in [";
-    for (unsigned int i=0;i<get_number_of_bins();++i) msg()<<get_bin_low_edge(i)<<", ";
-    msg()<<get_bin_up_edge(get_number_of_bins()-1)<<"]"<< endmsg;
+    ATH_MSG(INFO) << optprint <<"  "<<get_number_of_bins()<<" functions : ";
+    for (unsigned int i=0;i<get_number_of_bins();++i) msg()<<get_bin_low_edge(i)<<" < ("<<get_function(i)<<") < ";
+    msg()<<get_bin_up_edge(get_number_of_bins()-1)<< endmsg;
   }  
 }
 

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCSLateralShapeParametrizationHitChain.cxx

@@ -52,38 +52,47 @@ FCSReturnCode TFCSLateralShapeParametrizationHitChain::simulate(TFCSSimulationSt
     return FCSFatal;
   }
 
-  float Ehit=simulstate.E(calosample())/nhit;
+  float Elayer=simulstate.E(calosample());
+  float Ehit=Elayer/nhit;
+  float sumEhit=0;
 
   bool debug = msgLvl(MSG::DEBUG);
   if (debug) {
-    ATH_MSG_DEBUG("E("<<calosample()<<")="<<simulstate.E(calosample())<<" #hits="<<nhit);
+    ATH_MSG_DEBUG("E("<<calosample()<<")="<<simulstate.E(calosample())<<" #hits~"<<nhit);
   }
 
-  for (int i = 0; i < nhit; ++i) {
-    TFCSLateralShapeParametrizationHitBase::Hit hit; 
+  int ihit=0;
+  TFCSLateralShapeParametrizationHitBase::Hit hit; 
+  do {
+    hit.reset();
     hit.E()=Ehit;
     for(TFCSLateralShapeParametrizationHitBase* hitsim : m_chain) {
       if (debug) {
-        if (i < 2) hitsim->setLevel(MSG::DEBUG);
+        if (ihit < 2) hitsim->setLevel(MSG::DEBUG);
         else hitsim->setLevel(MSG::INFO);
       }
 
       for (int i = 0; i <= FCS_RETRY_COUNT; i++) {
+        //TODO: potentially change logic in case of a retry to redo the whole hit chain from an empty hit instead of just redoing one step in the hit chain
         if (i > 0) ATH_MSG_WARNING("TFCSLateralShapeParametrizationHitChain::simulate(): Retry simulate_hit call " << i << "/" << FCS_RETRY_COUNT);
   
         FCSReturnCode status = hitsim->simulate_hit(hit, simulstate, truth, extrapol);
 
-        if (status == FCSSuccess)
+        if (status == FCSSuccess) {
+          if(sumEhit+hit.E()>Elayer) hit.E()=Elayer-sumEhit;//sum of all hit energies needs to be Elayer: correct last hit accordingly
           break;
-        else if (status == FCSFatal)
-          return FCSFatal;
+        } else {
+          if (status == FCSFatal) return FCSFatal;
+        }    
 
         if (i == FCS_RETRY_COUNT) {
           ATH_MSG_ERROR("TFCSLateralShapeParametrizationHitChain::simulate(): simulate_hit call failed after " << FCS_RETRY_COUNT << "retries");
         }
       }
     }
-  }
+    sumEhit+=hit.E();
+    ++ihit;
+  } while (sumEhit<Elayer);
 
   return FCSSuccess;
 }