diff --git a/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSD.cxx b/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSD.cxx
index 2dbe767cab13b00f36d76d91660dfa36236edc52..b1d1dd2031d91136e332079d8c3df97a37a7f68d 100644
--- a/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSD.cxx
+++ b/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);
diff --git a/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSD.h b/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSD.h
index 3e4db72dec2759e8e9c7b67bdadb5fc280d0204f..d558c0e74e1a55a1e0b3680d7f4831dc3ba66341 100644
--- a/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSD.h
+++ b/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
diff --git a/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSDTool.cxx b/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSDTool.cxx
index 0a82737011ee04aa5412ba586ca536af5ece0de4..f2715e28a99a32d357948f4d9fb0dbbb1359ca30 100644
--- a/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSDTool.cxx
+++ b/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;
}
diff --git a/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSDTool.h b/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSDTool.h
index ceb5b94581a3a04b4b708b7d21cafc45ce3548a7..2ddc736d53fda5285c0310134b31050f925f96cd 100644
--- a/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSDTool.h
+++ b/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
diff --git a/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.cxx b/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.cxx
index d63b9c4c07c564d563559b77351868751bf12363..f99ce4b4a9d68c0df89d918aa96d26d53a96efbc 100644
--- a/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.cxx
+++ b/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.cxx
@@ -15,20 +15,20 @@ StatusCode CaloSimHitAlg::initialize()
// 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;Fraction;Events", 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_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_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_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;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;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));