diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt index bcbfe2e5a0322d4f85b14de8839386992e2e69b3..c4ee8f78b1ae3ddbd5873a0fdc81a645bb4897e8 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/CMakeLists.txt +++ b/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 diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/LinkDef.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/LinkDef.h index 6945bf8c0f9fc5a5eefc59e230859cf9c54fdcc6..53260840e54148fdfbc951de30eb62dae61d0e9d 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/LinkDef.h +++ b/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+; diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS2DFunctionLateralShapeParametrization.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS2DFunctionLateralShapeParametrization.h new file mode 100644 index 0000000000000000000000000000000000000000..f5f0af2b98cb71384745997976b815c1d6f4b292 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/ISF_FastCaloSimEvent/TFCS2DFunctionLateralShapeParametrization.h @@ -0,0 +1,66 @@ +/* + 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 diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCS2DFunctionLateralShapeParametrization.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCS2DFunctionLateralShapeParametrization.cxx new file mode 100644 index 0000000000000000000000000000000000000000..89b895f0201e74c4a744d0bb287cae812f781c12 --- /dev/null +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimEvent/src/TFCS2DFunctionLateralShapeParametrization.cxx @@ -0,0 +1,150 @@ +/* + 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)"); + } + } +}