Merge branch 'master-22.0.29-ecal' into 'master'

Master 22.0.29 ecal

See merge request faser/calypso!126

Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSD.cxx

@@ -52,6 +52,7 @@ G4bool EcalSensorSD::ProcessHits(G4Step* aStep, G4TouchableHistory* /*ROhist*/)
       return false;
   }
   edep *= CLHEP::MeV;
+
   //
   // Get the Touchable History:
   //
@@ -101,6 +102,23 @@ G4bool EcalSensorSD::ProcessHits(G4Step* aStep, G4TouchableHistory* /*ROhist*/)
   lP2[1] = localPosition2[1]*CLHEP::mm;
   lP2[0] = localPosition2[0]*CLHEP::mm;
 
+  // Birk's Law Correction
+
+  const G4Track*              const track    = aStep->GetTrack      () ;
+  const G4ParticleDefinition* const particle = track->GetDefinition () ;
+  //const double                      charge   = particle->GetPDGCharge  () ; 
+  const G4StepPoint* const          preStep  = aStep->GetPreStepPoint () ;
+  const G4MaterialCutsCouple* const material = preStep->GetMaterialCutsCouple () ;
+
+
+  //double ecorrected = edep;
+
+  double ecorrected = edep * birkCorrection( particle , track->GetKineticEnergy(), material ) ;
+
+  // Non-uniformity correction
+
+  ecorrected *= localNonUniformity(lP1[0], lP1[1]) ;
+
   // Now Navigate the history to know in what detector the step is:
   // and finally set the ID of det element in which the hit is.
   //
@@ -115,7 +133,7 @@ G4bool EcalSensorSD::ProcessHits(G4Step* aStep, G4TouchableHistory* /*ROhist*/)
   FaserTrackHelper trHelp(aStep->GetTrack());
   m_HitColl->Emplace(lP1,
                      lP2,
-                     edep,
+                     ecorrected,
                      aStep->GetPreStepPoint()->GetGlobalTime(),//use the global time. i.e. the time from the beginning of the event
                      trHelp.GetParticleLink(),
                      row,module);

Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSD.h

@@ -20,6 +20,11 @@
 class G4Step;
 class G4TouchableHistory;
 
+// G4 needed classes
+#include "Geant4/G4EnergyLossTables.hh"
+#include "Geant4/G4Material.hh"
+#include "Geant4/G4MaterialCutsCouple.hh"
+
 class EcalSensorSD : public G4VSensitiveDetector
 {
 public:
@@ -40,13 +45,268 @@ public:
       same SD classes as the standard simulation. */
   template <class... Args> void AddHit(Args&&... args){ m_HitColl->Emplace( args... ); }
 
+  /// Set the first coefficient of Birks's law   
+  void birk_c1 (double c1) { m_birk_c1 = c1 ; }   
+  /// Set the second coefficient of Birks's law   
+  void birk_c2 (double c2) { m_birk_c2 = c2; }   
+  /// Set the correction to the first coefficient of Birks's law   
+  void birk_c1cor (double c1cor) { m_birk_c1correction = c1cor; }
+
+  /// Amplitudes of the local non uniformity correction
+  void local_non_uniformity (double local) { m_a_local_outer_ecal = local ; }
+  /// Amplitudes of the global non uniformity correction
+  void global_non_uniformity (double global) { m_a_global_outer_ecal = global ; }
+
+  /// Correction for light reflection at the edges
+  void reflection_height(double h) { m_a_reflection_height = h ; }
+  void reflection_width(double w) { m_a_reflection_width = w ; }
+
 private:
   void indexMethod(const G4TouchableHistory *myTouch, int &station, int &plate);
   static const uint32_t kModuleDepth { 3 };
   static const uint32_t kFiberDepth  { 5 };
+
 protected:
   // The hits collection
   SG::WriteHandle<CaloHitCollection> m_HitColl;
+
+private:
+  // the first coefficient of Birk's law                    (c1)   
+  double m_birk_c1             ;   
+  // the second coefficient of Birk's law                   (c2)   
+  double m_birk_c2             ;   
+  // the correction to the first coefficient of Birk's law  (c1')   
+  double m_birk_c1correction   ;
+  
+  // Amplitudes of the local non uniformity correction
+  double m_a_local_outer_ecal  ;
+  // Amplitudes of the global non uniformity correction
+  double m_a_global_outer_ecal  ;
+  // Correction for light reflection at the edges
+  double  m_a_reflection_height ;
+  double  m_a_reflection_width  ;
+
+protected:
+  /// the first coefficient of Birks's law   
+  inline double birk_c1 () const { return m_birk_c1 ; }   
+  /// the second coefficient of Birks's law   
+  inline double birk_c2 () const { return m_birk_c2 ; }   
+  /// the correction to the first coefficient of Birks's law   
+  inline double birk_c1cor () const   { return m_birk_c1correction ; }
+  
+  /** Correction factor from Birk's Law
+   *  @param step current G4step
+   *  @return the correction factor
+   */
+  inline double birkCorrection ( const G4Step* step ) const ;
+  
+  /** Birk's correction for given particle with given kinetic energy
+   *  for the given material
+   *  @param  particle pointer to particle definition
+   *  @param  Ekine    particle kinetic energy
+   *  @param  maerial  pointer ot teh material
+   */
+  inline double  birkCorrection
+  ( const G4ParticleDefinition* particle ,
+    const double                eKine    ,
+    const G4MaterialCutsCouple* material ) const ;
+  
+  /** evaluate the correction for Birk's law
+   *  @param charge   the charge of the particle
+   *  @param dEdX     the nominal dEdX in the material
+   *  @param material the pointer ot teh material
+   *  @return the correction coefficient
+   */
+  inline double  birkCorrection
+  ( const double      charge   ,
+    const double      dEdX     ,
+    const G4Material* material ) const ;
+  
+  /** evaluate the correction for Birk's law
+   *  @param charge   the charge of the particle
+   *  @param dEdX     the nominal dEdX in the material
+   *  @param density  the density ot the material
+   *  @return the correction coefficient
+   */
+  inline double  birkCorrection
+  ( const double      charge   ,
+    const double      dEdX     ,
+    const double      density  ) const ;
+  
+  /** Correction due to the local non uniformity due to the light
+   *  collection efficiency in cell
+   */
+  inline double localNonUniformity (double x, double y) const ;
+
+  /** The fractions of energy deposited in consequitive time-bins
+   *  in the given Ecal/Hcal cell
+   *  @see CaloSensDet
+   *  @param time global time of energy deposition
+   *  @param cell cellID of the cell
+   *  @param slot (out) the first time bin
+   *  @param fracs (out) the vector of fractions for subsequent time-slots;
+   *  @return StatuscCode
+   */
+//   StatusCode timing
+//   ( const double            time      ,
+//     Fractions& fractions ) const;
+
+  
 };
 
+// ============================================================================
+/** evaluate the correction for Birk's law
+ *  @param charge   the charge of the particle
+ *  @param dEdX     the nominal dEdX in the material
+ *  @param density
+ *  @return the correction coefficient
+ */
+// ============================================================================
+inline double  EcalSensorSD::birkCorrection
+( const double      charge   ,
+  const double      dEdX     ,
+  const double      density  ) const
+{
+  const double C1 =
+    fabs( charge )  < 1.5 ? birk_c1() : birk_c1() * birk_c1cor() ;
+  const double C2 = birk_c2() ;
+
+  const double alpha = dEdX/ density ;
+
+  return 1.0 / ( 1.0 + C1 * alpha + C2 * alpha * alpha ) ;
+}
+
+// ============================================================================
+/** evaluate the correction for Birk's law
+ *  @param charge   the charge of the particle
+ *  @param dEdX     the nominal dEdX in the material
+ *  @param material the pointer ot teh material
+ *  @return the correction coefficient
+ */
+// ============================================================================
+inline double  EcalSensorSD::birkCorrection
+( const double      charge   ,
+  const double      dEdX     ,
+  const G4Material* material ) const
+{
+  if ( 0 == material )
+    { std::cerr << "birkCorrection(): invalid material " << std::endl; return 1. ; } // RETURN
+  return birkCorrection( charge , dEdX , material->GetDensity() ) ;
+}
+
+// ============================================================================
+/** Birk's correction for given particle with given kinetic energy
+ *  for the given material
+ *  @param  particle pointer to particle definition
+ *  @param  Ekine    particle kinetic energy
+ *  @param  maerial  pointer ot teh material
+ */
+// ============================================================================
+inline double EcalSensorSD::birkCorrection
+( const G4ParticleDefinition* particle ,
+  const double                eKine    ,
+  const G4MaterialCutsCouple* material ) const
+{
+  if (  0 == particle || 0 == material )
+    { std::cerr << "birkCorrection(): invalid parameters " << std::endl; return 1.0 ; } // RETURN
+
+  const double charge = particle -> GetPDGCharge() ;
+  if ( fabs ( charge ) < 0.1                      ) { return 1.0 ; } // EEUTRN
+
+  // get the nominal dEdX
+  const double dEdX  =
+    G4EnergyLossTables::GetDEDX ( particle , eKine , material ) ;
+
+  // make an actual evaluation
+  return birkCorrection
+    ( charge ,
+      dEdX   ,
+      material->GetMaterial()->GetDensity() ) ;
+}
+
+// ============================================================================
+// Birk's Law
+// ============================================================================
+/** Correction factor from Birk's Law
+ *  Factor = 1/(1+C1*dEdx/rho+C2*(dEdx/rho)^2)
+ *  Where :
+ *      - C1 = 0.013 g/MeV/cm^2 [Ref NIM 80 (1970) 239]
+ *      - C2 = 9.6.10^-6 g^2/MeV^2/cm^4
+ *      - dEdx in MeV/cm
+ *      - rho = density in g/cm^3
+ */
+// ============================================================================
+inline double EcalSensorSD::birkCorrection ( const G4Step* step ) const
+{
+  if ( !step ) { return 1 ; }                               // RETURN
+
+  // Track
+  const G4Track* track  = step->GetTrack() ;
+  const double charge   = track->GetDefinition()->GetPDGCharge()  ;
+
+  // Only for charged tracks
+  if ( fabs( charge ) < 0.1 ) { return 1 ; }                 // RETURN
+
+  // make an actual evaluation
+  return birkCorrection
+    ( track->GetDefinition         () ,
+      track->GetKineticEnergy      () ,
+      track->GetMaterialCutsCouple () ) ;
+}
+// ============================================================================
+
+inline double EcalSensorSD::localNonUniformity (double x, double y) const
+{
+
+  // Only for ECal for the moment
+  double correction = 1. ;
+
+  // Find the position of the step relative to the centre of the cell
+  // Done before passing in so set centres x0 = y0 = 0
+  double x0       = 0;
+  double y0       = 0;
+
+
+  // Distance between fibers
+  // and correction amplitude
+  double d        = 15.25 * CLHEP::mm ;
+  double A_local  = m_a_local_outer_ecal ; 
+  double A_global = m_a_global_outer_ecal ;
+
+  // Cell size
+  double cellSize = 121.9 * CLHEP::mm; 
+
+  float lowTolerance = 1e-20;
+
+  // Local uniformity is product of x and y sine-like functions
+  // The Amplitude of the sin-like function is a function of x and y
+  if ( A_local > lowTolerance )
+    correction += A_local / 2. * ( 1. - cos( 2. * CLHEP::pi * (x-x0)/d ) ) *
+      ( 1. - cos( 2. * CLHEP::pi * (y-y0)/d ) ) ;
+
+  double rX(0.) , rY(0.) , hCell(0.) ;
+
+  // Global non uniformity
+  if ( A_global > lowTolerance ) {
+    rX = x - x0 ;
+    rY = y - y0 ;
+    hCell = cellSize / 2. ;
+    correction +=
+      A_global * ( hCell - rX ) * ( rX + hCell ) / ( hCell * hCell )
+      * ( hCell - rY ) * ( rY + hCell ) / ( hCell * hCell ) ;
+  }
+
+  // Light Reflexion on the edges
+  if ( m_a_reflection_height > lowTolerance ) {
+    rX = rX / m_a_reflection_width ;
+    rY = rY / m_a_reflection_width ;
+    hCell = hCell / m_a_reflection_width ;
+    correction += m_a_reflection_height *
+      ( exp( - fabs ( rX + hCell ) ) + exp( - fabs ( rX - hCell ) )
+        + exp( - fabs ( rY + hCell ) ) + exp( - fabs ( rY - hCell ) ) ) ;
+  }
+
+  return correction ;
+}
+
 #endif //ECALG4_SD_ECALSENSORSD_H

Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSDTool.cxx

@@ -18,8 +18,23 @@
 
 EcalSensorSDTool::EcalSensorSDTool(const std::string& type, const std::string& name, const IInterface* parent)
   : SensitiveDetectorBase( type , name , parent )
+  , m_birk_c1 ( 0.013  * CLHEP::g/CLHEP::MeV/CLHEP::cm2 ) // 1st coef. of Birk's   
+  , m_birk_c2 ( 9.6E-6 * CLHEP::g*CLHEP::g/CLHEP::MeV/CLHEP::MeV/CLHEP::cm2/CLHEP::cm2 ) // 2nd coef. of Birk's law 
+  , m_birk_c1correction  ( 0.57142857)  //correction of c1 for 2e charged part.  
+  , m_a_local_outer_ecal ( 0. )
+  , m_a_global_outer_ecal  ( 0.03  ) // global non uniformity amplitude
+  , m_a_reflection_height ( 0.09 ) // reflection on the edges - height
+  , m_a_reflection_width  ( 6. * CLHEP::mm ) // reflection on the edges - width
 {
 
+  declareProperty ( "BirkC1"               ,  m_birk_c1               ) ;   
+  declareProperty ( "BirkC1cor"            ,  m_birk_c1correction     ) ;   
+  declareProperty ( "BirkC2"               ,  m_birk_c2               ) ;
+  declareProperty ( "LocalNonUnifomity"    ,  m_a_local_outer_ecal  ) ;
+  declareProperty ( "GlobalNonUnifomity"   ,  m_a_global_outer_ecal ) ;
+  declareProperty ( "ReflectionHeight"     ,  m_a_reflection_height   ) ;  
+  declareProperty ( "ReflectionWidth"      ,  m_a_reflection_width    ) ;
+
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -27,6 +42,17 @@ EcalSensorSDTool::EcalSensorSDTool(const std::string& type, const std::string& n
 G4VSensitiveDetector* EcalSensorSDTool::makeSD() const
 {
   ATH_MSG_DEBUG( "Creating Ecal SD: " << name() );
-  return new EcalSensorSD(name(), m_outputCollectionNames[0]);
+
+  EcalSensorSD* ecsd = new EcalSensorSD(name(), m_outputCollectionNames[0]);
+
+  ecsd->birk_c1(m_birk_c1);
+  ecsd->birk_c1cor(m_birk_c1correction);
+  ecsd->birk_c2(m_birk_c2);
+  ecsd->local_non_uniformity(m_a_local_outer_ecal);
+  ecsd->global_non_uniformity(m_a_global_outer_ecal);
+  ecsd->reflection_height(m_a_reflection_height);
+  ecsd->reflection_width(m_a_reflection_width);
+
+  return ecsd;
 }
 

Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSDTool.h

@@ -21,7 +21,8 @@ class G4VSensitiveDetector;
 
 class EcalSensorSDTool : public SensitiveDetectorBase
 {
- public:
+
+public:
   // Constructor
   EcalSensorSDTool(const std::string& type, const std::string& name, const IInterface *parent);
 
@@ -31,6 +32,22 @@ class EcalSensorSDTool : public SensitiveDetectorBase
 protected:
   // Make me an SD!
   G4VSensitiveDetector* makeSD() const override final;
+
+private:
+  // the first coefficient of Birk's law                    (c1)   
+  double m_birk_c1             ;   
+  // the second coefficient of Birk's law                   (c2)   
+  double m_birk_c2             ;   
+  // the correction to the first coefficient of Birk's law  (c1')   
+  double m_birk_c1correction   ;
+  
+  // Amplitudes of the local non uniformity correction
+  double m_a_local_outer_ecal  ;
+  // Amplitudes of the global non uniformity correction
+  double m_a_global_outer_ecal  ;
+  // Correction for light reflection at the edges
+  double  m_a_reflection_height ;
+  double  m_a_reflection_width  ;
 };
 
 #endif //ECALG4_SD_ECALSENSORSDTOOL_H

Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.cxx

@@ -9,26 +9,26 @@ CaloSimHitAlg::~CaloSimHitAlg() { }
 
 StatusCode CaloSimHitAlg::initialize()
 
-{
+{ 
 
     // initialize a histogram 
     // letter at end of TH1 indicated variable type (D double, F float etc)
     // first string is root object name, second is histogram title
 
-    m_eloss = new TH1D("eLoss", "Calo Hit Energy Loss fraction", 1000, 0, 0.0001); 
-    m_elossTL = new TH1D("eLossTL", "Calo Hit Energy Loss Fraction Top Left Module", 1000, 0, 0.0001);
-    m_elossTR = new TH1D("eLossTR", "Calo Hit Energy Loss Fraction Top Right Module", 1000, 0, 0.0001);
-    m_elossBR = new TH1D("eLossBR", "Calo Hit Energy Loss Fraction Bottom Right Module", 1000, 0, 0.0001);
-    m_elossBL = new TH1D("eLossBL", "Calo Hit Energy Loss Fraction Bottom Left Module", 1000, 0, 0.0001);
+    m_eloss = new TH1D("eLoss", "Calo Hit Energy Loss Fraction;Fraction;Events", 1000, 0, 0.0001); 
+    m_elossTL = new TH1D("eLossTL", "Calo Hit Energy Loss Fraction Top Left Module;Fraction;Events", 1000, 0, 0.0001);
+    m_elossTR = new TH1D("eLossTR", "Calo Hit Energy Loss Fraction Top Right Module;Fraction;Events", 1000, 0, 0.0001);
+    m_elossBR = new TH1D("eLossBR", "Calo Hit Energy Loss Fraction Bottom Right Module;Fraction;Events", 1000, 0, 0.0001);
+    m_elossBL = new TH1D("eLossBL", "Calo Hit Energy Loss Fraction Bottom Left Module;Fraction;Events", 1000, 0, 0.0001);
 
-    m_elossTot = new TH1D("eLossTot", "Total Energy Fraction Deposited in Calorimeter", 1000, 0, 1);   
-    m_elossTLTot = new TH1D("eLossTLTot", "Total Energy Fraction Deposited in Top Left Module", 1000, 0, 1);
-    m_elossTRTot = new TH1D("eLossTRTot", "Total Energy Fraction Deposited in Top Right Module", 1000, 0, 1);
-    m_elossBRTot = new TH1D("eLossBRTot", "Total Energy Fraction Deposited in Bottom Right Module", 1000, 0, 1);
-    m_elossBLTot = new TH1D("eLossBLTot", "Total Energy Fraction Deposited in Bottom Left Module", 1000, 0, 1);    
+    m_elossTot = new TH1D("eLossTot", "Total Energy Fraction Deposited in Calorimeter;Fraction;Events", 1000, 0, 1);   
+    m_elossTLTot = new TH1D("eLossTLTot", "Total Energy Fraction Deposited in Top Left Module;Fraction;Events", 1000, 0, 1);
+    m_elossTRTot = new TH1D("eLossTRTot", "Total Energy Fraction Deposited in Top Right Module;Fraction;Events", 1000, 0, 1);
+    m_elossBRTot = new TH1D("eLossBRTot", "Total Energy Fraction Deposited in Bottom Right Module;Fraction;Events", 1000, 0, 1);
+    m_elossBLTot = new TH1D("eLossBLTot", "Total Energy Fraction Deposited in Bottom Left Module;Fraction;Events", 1000, 0, 1);    
 
-    m_module = new TH2D("module", "Calo Hit Module", 2, -0.5, 1.5, 2, -0.5, 1.5 );
-    m_modulePos = new TH2D("modulePos", "Calo Hit Module Position", 13, -130, 130, 13, -130, 130 );
+    m_module = new TH2D("module;x;y", "Calo Hit Module", 2, -0.5, 1.5, 2, -0.5, 1.5 );
+    m_modulePos = new TH2D("modulePos;x;y", "Calo Hit Module Position", 13, -130, 130, 13, -130, 130 );
         
 
     ATH_CHECK(histSvc()->regHist("/HIST/eloss", m_eloss));