Add 2D lateral shape simulation based on general 2D functions

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt

@@ -68,6 +68,7 @@ atlas_add_root_dictionary( ISF_FastCaloSimEvent _dictSource
                            ISF_FastCaloSimEvent/TFCSFlatLateralShapeParametrization.h
                            ISF_FastCaloSimEvent/TFCSHistoLateralShapeParametrization.h
                            ISF_FastCaloSimEvent/TFCSHistoLateralShapeParametrizationFCal.h
+                           ISF_FastCaloSimEvent/TFCS2DFunctionLateralShapeParametrization.h
                            ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h
                            ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeightHitAndMiss.h
                            ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitNumberFromE.h

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

@@ -47,6 +47,7 @@
 #include "ISF_FastCaloSimEvent/TFCSFlatLateralShapeParametrization.h"
 #include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeParametrization.h"
 #include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeParametrizationFCal.h"
+#include "ISF_FastCaloSimEvent/TFCS2DFunctionLateralShapeParametrization.h"
 #include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h"
 #include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeightHitAndMiss.h"
 #include "ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitNumberFromE.h"
@@ -205,6 +206,7 @@
 #pragma link C++ class TFCSFlatLateralShapeParametrization+;
 #pragma link C++ class TFCSHistoLateralShapeParametrization+;
 #pragma link C++ class TFCSHistoLateralShapeParametrizationFCal+;
+#pragma link C++ class TFCS2DFunctionLateralShapeParametrization+;
 #pragma link C++ class TFCSHistoLateralShapeWeight+;
 #pragma link C++ class TFCSHistoLateralShapeWeightHitAndMiss+;
 #pragma link C++ class TFCSLateralShapeParametrizationHitNumberFromE+;

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS2DFunctionLateralShapeParametrization.h

+/*
+  Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef TFCS2DFunctionLateralShapeParametrization_h
+#define TFCS2DFunctionLateralShapeParametrization_h
+
+#include "ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitBase.h"
+#include "ISF_FastCaloSimEvent/TFCS2DFunction.h"
+#include "ISF_FastCaloSimEvent/TFCSTruthState.h"
+
+class TH2;
+
+class TFCS2DFunctionLateralShapeParametrization:public TFCSLateralShapeParametrizationHitBase {
+public:
+  TFCS2DFunctionLateralShapeParametrization(const char* name=nullptr, const char* title=nullptr);
+  ~TFCS2DFunctionLateralShapeParametrization();
+
+  ///Status bit for FCS needs
+  enum FCSStatusBits {
+     k_phi_symmetric = BIT(15) ///< Set this bit to simulate phi symmetric histograms
+  };
+
+  bool is_phi_symmetric() const {return TestBit(k_phi_symmetric);};
+  virtual void set_phi_symmetric() {SetBit(k_phi_symmetric);};
+  virtual void reset_phi_symmetric() {ResetBit(k_phi_symmetric);};
+
+  /// set the desired number of hits
+  void set_number_of_hits(float nhits);
+
+  float get_number_of_expected_hits() const {return m_nhits;};
+  
+  ///default for this class is to simulate get_number_of_expected_hits() hits, 
+  ///which gives fluctuations sigma^2=1/get_number_of_expected_hits()
+  virtual double get_sigma2_fluctuation(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const override;
+
+  /// default for this class is to simulate get_number_of_expected_hits() hits
+  int get_number_of_hits(TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) const override;
+
+  /// 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 FCSReturnCode simulate_hit(Hit& hit,TFCSSimulationState& simulstate,const TFCSTruthState* truth, const TFCSExtrapolationState* extrapol) override;
+
+  /// Init from function
+  bool Initialize(TFCS2DFunction* func,float nhits=-1);
+  
+  TFCS2DFunction* function() {return m_function;};
+  const TFCS2DFunction* function() const {return m_function;};
+  
+  void Print(Option_t *option = "") const override;
+protected:
+  /// Histogram to be used for the shape simulation
+  TFCS2DFunction* m_function;
+  float m_nhits;
+
+private:
+
+  ClassDefOverride(TFCS2DFunctionLateralShapeParametrization,2)  //TFCS2DFunctionLateralShapeParametrization
+};
+
+#if defined(__ROOTCLING__) && defined(__FastCaloSimStandAlone__)
+#pragma link C++ class TFCS2DFunctionLateralShapeParametrization+;
+#endif
+
+#endif

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

+/*
+  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/TFCS2DFunctionLateralShapeParametrization.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 "TH2.h"
+
+#include "HepPDT/ParticleData.hh"
+#include "HepPDT/ParticleDataTable.hh"
+
+//=============================================
+//======= TFCS2DFunctionLateralShapeParametrization =========
+//=============================================
+
+TFCS2DFunctionLateralShapeParametrization::TFCS2DFunctionLateralShapeParametrization(const char* name, const char* title) :
+  TFCSLateralShapeParametrizationHitBase(name,title),m_function(nullptr),m_nhits(0)
+{
+  reset_phi_symmetric();
+}
+
+TFCS2DFunctionLateralShapeParametrization::~TFCS2DFunctionLateralShapeParametrization()
+{
+  if(m_function) delete m_function;
+  m_function=nullptr;
+}
+
+double TFCS2DFunctionLateralShapeParametrization::get_sigma2_fluctuation(TFCSSimulationState& /*simulstate*/,const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* /*extrapol*/) const
+{
+  return 1.0/m_nhits;
+}
+
+int TFCS2DFunctionLateralShapeParametrization::get_number_of_hits(TFCSSimulationState &simulstate, const TFCSTruthState* /*truth*/, const TFCSExtrapolationState* /*extrapol*/) const
+{
+  if (!simulstate.randomEngine()) {
+    return -1;
+  }
+
+  return CLHEP::RandPoisson::shoot(simulstate.randomEngine(), m_nhits);
+}
+
+void TFCS2DFunctionLateralShapeParametrization::set_number_of_hits(float nhits)
+{
+  m_nhits=nhits;
+}
+
+FCSReturnCode TFCS2DFunctionLateralShapeParametrization::simulate_hit(Hit &hit, TFCSSimulationState &simulstate, const TFCSTruthState* truth, const TFCSExtrapolationState* /*extrapol*/)
+{
+  if (!simulstate.randomEngine()) {
+    return FCSFatal;
+  }
+  if (!m_function) {
+    return FCSFatal;
+  }
+  
+  const int     pdgId    = truth->pdgid();
+  const double  charge   = HepPDT::ParticleID(pdgId).charge();
+
+  const int cs=calosample();
+  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();
+
+  if (TMath::IsNaN(center_r) or TMath::IsNaN(center_z) or TMath::IsNaN(center_eta) or TMath::IsNaN(center_phi)) { //Check if extrapolation fails
+    return FCSFatal;
+  }
+
+  float alpha, r, rnd1, rnd2;
+  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;
+      rnd2*=2;
+      m_function->rnd_to_fct(alpha,r,rnd1,rnd2);
+      alpha=-alpha;
+    } else { //Fill positive phi half of shape
+      rnd2*=2;
+      m_function->rnd_to_fct(alpha,r,rnd1,rnd2);
+    }
+  } else {
+    m_function->rnd_to_fct(alpha,r,rnd1,rnd2);
+  }
+  if(TMath::IsNaN(alpha) || TMath::IsNaN(r)) {
+    ATH_MSG_ERROR("  2D function, #hits="<<m_nhits<<" alpha="<<alpha<<" r="<<r<<" rnd1="<<rnd1<<" rnd2="<<rnd2);
+    alpha=0;
+    r=0.001;
+
+    ATH_MSG_ERROR("  This error could probably be retried");
+    return FCSFatal;
+  }
+  
+  float delta_eta_mm = r * cos(alpha);
+  float delta_phi_mm = r * sin(alpha);
+  
+  // Particles with negative eta are expected to have the same shape as those with positive eta after transformation: delta_eta --> -delta_eta
+  if(center_eta<0.)delta_eta_mm = -delta_eta_mm;
+  // Particle with negative charge are expected to have the same shape as positively charged particles after transformation: delta_phi --> -delta_phi
+  if(charge < 0.)  delta_phi_mm = -delta_phi_mm;
+
+  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 = delta_eta_mm / eta_jakobi / dist000;
+  const float delta_phi = delta_phi_mm / center_r;
+
+  hit.setEtaPhiZE(center_eta + delta_eta,center_phi + delta_phi,center_z, hit.E());
+
+  ATH_MSG_DEBUG("HIT: E="<<hit.E()<<" cs="<<cs<<" eta="<<hit.eta()<<" phi="<<hit.phi()<< " z="<<hit.z()<<" r="<<r<<" alpha="<<alpha);
+
+  return FCSSuccess;
+}
+
+
+bool TFCS2DFunctionLateralShapeParametrization::Initialize(TFCS2DFunction* func,float nhits)
+{
+  if(!func) return false;
+  if(m_function) delete m_function;
+  m_function=func;
+
+  if(nhits>0) set_number_of_hits(nhits);
+
+  return true;
+}
+
+void TFCS2DFunctionLateralShapeParametrization::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(is_phi_symmetric()) {
+      ATH_MSG_INFO(optprint <<"  2D function, #hits="<<m_nhits<<" (phi symmetric)");
+    } else {
+      ATH_MSG_INFO(optprint <<"  2D function, #hits="<<m_nhits<<" (not phi symmetric)");
+    }
+  }  
+}