add functionality for Gaus(Log(E)) weighted hits

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt

@@ -74,7 +74,9 @@ atlas_add_root_dictionary( ISF_FastCaloSimEvent _dictSource
                            ISF_FastCaloSimEvent/TFCSHistoLateralShapeParametrizationFCal.h
                            ISF_FastCaloSimEvent/TFCS2DFunctionLateralShapeParametrization.h
                            ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h
+                           ISF_FastCaloSimEvent/TFCSHistoLateralShapeGausLogWeight.h
                            ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeightHitAndMiss.h
+                           ISF_FastCaloSimEvent/TFCSHistoLateralShapeGausLogWeightHitAndMiss.h
                            ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitNumberFromE.h
                            ISF_FastCaloSimEvent/TFCSHitCellMapping.h
                            ISF_FastCaloSimEvent/TFCSHitCellMappingFCal.h

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

@@ -50,7 +50,9 @@
 #include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeParametrizationFCal.h"
 #include "ISF_FastCaloSimEvent/TFCS2DFunctionLateralShapeParametrization.h"
 #include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h"
+#include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeGausLogWeight.h"
 #include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeightHitAndMiss.h"
+#include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeGausLogWeightHitAndMiss.h"
 #include "ISF_FastCaloSimEvent/TFCSLateralShapeParametrizationHitNumberFromE.h"
 #include "ISF_FastCaloSimEvent/TFCSHitCellMapping.h"
 #include "ISF_FastCaloSimEvent/TFCSHitCellMappingFCal.h"
@@ -211,7 +213,9 @@
 #pragma link C++ class TFCSHistoLateralShapeParametrizationFCal+;
 #pragma link C++ class TFCS2DFunctionLateralShapeParametrization+;
 #pragma link C++ class TFCSHistoLateralShapeWeight+;
+#pragma link C++ class TFCSHistoLateralShapeGausLogWeight+;
 #pragma link C++ class TFCSHistoLateralShapeWeightHitAndMiss+;
+#pragma link C++ class TFCSHistoLateralShapeGausLogWeightHitAndMiss+;
 #pragma link C++ class TFCSLateralShapeParametrizationHitNumberFromE+;
 #pragma link C++ class TFCSHitCellMapping+;
 #pragma link C++ class TFCSHitCellMappingFCal+;

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSHistoLateralShapeGausLogWeight.h

+/*
+  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef TFCSHistoLateralShapeGausLogWeight_h
+#define TFCSHistoLateralShapeGausLogWeight_h
+
+#include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h"
+
+class TH1;
+
+class TFCSHistoLateralShapeGausLogWeight:public TFCSHistoLateralShapeWeight {
+public:
+  TFCSHistoLateralShapeGausLogWeight(const char* name=nullptr, const char* title=nullptr);
+  virtual ~TFCSHistoLateralShapeGausLogWeight();
+
+  /// 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) override;
+protected:
+
+  ClassDefOverride(TFCSHistoLateralShapeGausLogWeight,1)  //TFCSHistoLateralShapeGausLogWeight
+};
+
+#if defined(__ROOTCLING__) && defined(__FastCaloSimStandAlone__)
+#pragma link C++ class TFCSHistoLateralShapeGausLogWeight+;
+#endif
+
+#endif

Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCSHistoLateralShapeGausLogWeightHitAndMiss.h

+/*
+  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef TFCSHistoLateralShapeGausLogWeightHitAndMiss_h
+#define TFCSHistoLateralShapeGausLogWeightHitAndMiss_h
+
+#include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h"
+
+class TH1;
+
+class TFCSHistoLateralShapeGausLogWeightHitAndMiss:public TFCSHistoLateralShapeWeight {
+public:
+  TFCSHistoLateralShapeGausLogWeightHitAndMiss(const char* name=nullptr, const char* title=nullptr);
+  virtual ~TFCSHistoLateralShapeGausLogWeightHitAndMiss();
+
+  /// 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) override;
+
+private:
+
+  ClassDefOverride(TFCSHistoLateralShapeGausLogWeightHitAndMiss,1)  //TFCSHistoLateralShapeGausLogWeightHitAndMiss
+};
+
+#if defined(__ROOTCLING__) && defined(__FastCaloSimStandAlone__)
+#pragma link C++ class TFCSHistoLateralShapeGausLogWeightHitAndMiss+;
+#endif
+
+#endif

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

+/*
+  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "CLHEP/Random/RandFlat.h"
+#include "CLHEP/Random/RandGauss.h"
+
+#include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeGausLogWeight.h"
+#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
+
+#include "TH1.h"
+#include "TVector2.h"
+#include "TMath.h"
+
+//=============================================
+//======= TFCSHistoLateralShapeGausLogWeight =========
+//=============================================
+
+TFCSHistoLateralShapeGausLogWeight::TFCSHistoLateralShapeGausLogWeight(const char* name, const char* title) :
+  TFCSHistoLateralShapeWeight(name,title)
+{
+}
+
+TFCSHistoLateralShapeGausLogWeight::~TFCSHistoLateralShapeGausLogWeight()
+{
+}
+
+FCSReturnCode TFCSHistoLateralShapeGausLogWeight::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);
+  float RMS   =m_hist->GetBinError(bin);
+  if(RMS>0) {
+    float logweight=CLHEP::RandGauss::shoot(simulstate.randomEngine(), -0.5*RMS*RMS, RMS);
+    weight*=TMath::Exp(logweight);
+  }  
+  hit.E()*=weight;
+
+  ATH_MSG_DEBUG("HIT: E="<<hit.E()<<" dR_mm="<<delta_r_mm<<" weight="<<weight);
+  return FCSSuccess;
+}
+

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

+/*
+  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "CLHEP/Random/RandFlat.h"
+#include "CLHEP/Random/RandGauss.h"
+
+#include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeGausLogWeightHitAndMiss.h"
+#include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
+
+#include "TH1.h"
+#include "TVector2.h"
+#include "TMath.h"
+
+//=============================================
+//======= TFCSHistoLateralShapeGausLogWeightHitAndMiss =========
+//=============================================
+
+TFCSHistoLateralShapeGausLogWeightHitAndMiss::TFCSHistoLateralShapeGausLogWeightHitAndMiss(const char* name, const char* title):TFCSHistoLateralShapeWeight(name,title)
+{
+}
+
+TFCSHistoLateralShapeGausLogWeightHitAndMiss::~TFCSHistoLateralShapeGausLogWeightHitAndMiss()
+{
+}
+
+FCSReturnCode TFCSHistoLateralShapeGausLogWeightHitAndMiss::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 meanweight=m_hist->GetBinContent(bin);
+  float weight=meanweight;
+  float RMS   =m_hist->GetBinError(bin);
+  if(RMS>0) {
+    float logweight=CLHEP::RandGauss::shoot(simulstate.randomEngine(), -0.5*RMS*RMS, RMS);
+    weight*=TMath::Exp(logweight);
+  }  
+  if(meanweight<=1) {
+    //if meanweight<=1, give lower energy to hit. 
+    //TFCSLateralShapeParametrizationHitChain needs to be able to generate more hits in this case
+    hit.E()*=weight;
+  } else {
+    //if meanweight>1, accept only 1/meanweight events, but give them a higher energy increased by weight. 
+    //This leads to larger fluctuations, while keeping the shape unchanged.
+    float prob=1.0/meanweight;
+    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<<" meanweight="<<meanweight<<" weight="<<weight);
+  return FCSSuccess;
+}
+
+

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

@@ -3,6 +3,7 @@
 */
 
 #include "CLHEP/Random/RandFlat.h"
+#include "CLHEP/Random/RandGauss.h"
 
 #include "ISF_FastCaloSimEvent/TFCSHistoLateralShapeWeight.h"
 #include "ISF_FastCaloSimEvent/TFCSSimulationState.h"
@@ -51,6 +52,10 @@ FCSReturnCode TFCSHistoLateralShapeWeight::simulate_hit(Hit& hit,TFCSSimulationS
   if(bin<1) bin=1;
   if(bin>m_hist->GetNbinsX()) bin=m_hist->GetNbinsX();
   float weight=m_hist->GetBinContent(bin);
+  float RMS   =m_hist->GetBinError(bin);
+  if(RMS>0) {
+    weight=CLHEP::RandGauss::shoot(simulstate.randomEngine(), weight, RMS);
+  }  
   hit.E()*=weight;
 
   ATH_MSG_DEBUG("HIT: E="<<hit.E()<<" dR_mm="<<delta_r_mm<<" weight="<<weight);

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

@@ -49,31 +49,26 @@ FCSReturnCode TFCSHistoLateralShapeWeightHitAndMiss::simulate_hit(Hit& hit,TFCSS
   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);
+  float meanweight=m_hist->GetBinContent(bin);
+  float weight=meanweight;
   float RMS   =m_hist->GetBinError(bin);
   if(RMS>0) {
-    if(weight>=1) {
-      weight=CLHEP::RandGauss::shoot(simulstate.randomEngine(), weight, RMS);
-      //weight<1 can't be corrected for with HitAndMiss. So protect against a fluctuation leading to a weight<1
-      if(weight<1) weight=1;
-    } else {
-      weight=CLHEP::RandGauss::shoot(simulstate.randomEngine(), weight, RMS);
-    }
+    weight=CLHEP::RandGauss::shoot(simulstate.randomEngine(), meanweight, RMS);
   }  
-  if(weight<=1) {
-    //if weight<=1, give lower energy to hit. 
+  if(meanweight<=1) {
+    //if meanweight<=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. 
+    //if meanweight>1, accept only 1/meanweight 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 prob=1.0/meanweight;
     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);
+  ATH_MSG_DEBUG("HIT: E="<<hit.E()<<" dR_mm="<<delta_r_mm<<" meanweight="<<meanweight<<" weight="<<weight);
   return FCSSuccess;
 }