From d5619629280648b2db2baa2dc0aaa01175fdef71 Mon Sep 17 00:00:00 2001
From: Nicholas Styles <nicholas.styles@desy.de>
Date: Tue, 10 Nov 2020 11:44:00 +0000
Subject: [PATCH] Merge branch 'speedupOfMMStripsResponseSim' into '21.9'
Improve the CPU usage of the MM strips response simulation
See merge request atlas/athena!38002
---
.../MM_StripsResponseSimulation.h | 3 +-
.../src/MM_StripsResponseSimulation.cxx | 85 ++++++++++++-------
2 files changed, 55 insertions(+), 33 deletions(-)
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripsResponseSimulation.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripsResponseSimulation.h
index 2ed48255b45..e6deffa1d1a 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripsResponseSimulation.h
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripsResponseSimulation.h
@@ -182,7 +182,8 @@ private:
mutable Athena::MsgStreamMember m_msg = Athena::MsgStreamMember("MMStripResponseSimulation");
// seperate random number generation for performance monitoring
- float getTransversDiffusion(float posY);
+ float generateTransverseDiffusion(float posY);
+ float getTransverseDiffusion(float posY);
float getLongitudinalDiffusion(float posY);
float getEffectiveCharge();
float getPathLengthTraveled();
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripsResponseSimulation.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripsResponseSimulation.cxx
index 18216e7d3c9..7cfdbd7de16 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripsResponseSimulation.cxx
+++ b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripsResponseSimulation.cxx
@@ -233,23 +233,9 @@ void MM_StripsResponseSimulation::whichStrips( const float & hitx,
<< pathLengthTraveled
);
- for (auto& Electron : IonizationCluster.getElectrons()){
-
- // m_longitudinalDiffusionFunction->SetParameters(1.0, 0., initialPosition.Y()*m_longitudinalDiffusionSigma);
-
- // if ( m_longitudinalDiffusionSigma == 0 || m_transverseDiffusionSigma == 0) {
-
- // m_transverseDiffusionFunction->SetParameters( initialPosition.Y()*m_transverseDiffusionSigma , 0.0);
-
- // } else {
-
- // m_transverseDiffusionFunction->SetParameters( initialPosition.Y()*m_transverseDiffusionSigma , 1.0);
-
- // }
-
- Electron->setOffsetPosition(getTransversDiffusion(initialPosition.Y()) , getLongitudinalDiffusion(initialPosition.Y()) );
-
- }
+ for (auto& Electron : IonizationCluster.getElectrons()) {
+ Electron->setOffsetPosition(getTransverseDiffusion(initialPosition.Y()) , getLongitudinalDiffusion(initialPosition.Y()) );
+ }
IonizationCluster.propagateElectrons( lorentzAngle , m_driftVelocity );
@@ -337,28 +323,63 @@ void MM_StripsResponseSimulation::whichStrips( const float & hitx,
} // end of whichStrips()
-float MM_StripsResponseSimulation::getTransversDiffusion(float posY){
- if ( m_longitudinalDiffusionSigma == 0 || m_transverseDiffusionSigma == 0) {
- m_transverseDiffusionFunction->SetParameter(2, posY*m_transverseDiffusionSigma);
- m_transverseDiffusionFunction->SetParameter(5, 0.0);
- } else {
- m_transverseDiffusionFunction->SetParameter(2, posY*m_transverseDiffusionSigma);
- m_transverseDiffusionFunction->SetParameter(5, 1.0);
- }
- return m_transverseDiffusionFunction->GetRandom();
+float MM_StripsResponseSimulation::generateTransverseDiffusion(float posY) {
+ // this is a helper function used in getTransverseDiffusion, generating double gaussian distributions
+
+ // avoid division by zero in calculation of scale if ypos is 0
+ // also protect against negative values of ypos which should never happen
+ if (posY <= 0) posY = 0.001;
+
+ // need to scale weigths since initial distributions were not normalized
+ double scale = 0.001/(posY*m_transverseDiffusionSigma);
+
+ double uni = m_random->Uniform(0, 1.0+scale);
+ if (uni < scale) return m_random->Gaus(0.0, 1.0);
+ return m_random->Gaus(0., m_transverseDiffusionSigma * posY);
+}
+
+
+float MM_StripsResponseSimulation::getTransverseDiffusion(float posY) {
+ // the random numbers are generate from the following function:
+ // "1.*TMath::Exp(-TMath::Power(x,2.)/(2.*[0]*[0])) + 0.001*TMath::Exp(-TMath::Power(x,2)/(2.*[1]*[1]))"
+ // in the range from -1 to 1
+ // Since sampling from a TF1 where its parameters are changed for every random number is slow,
+ // the following approach is used to generate random numbers from the mixture of two Gaussian:
+ // https://stats.stackexchange.com/questions/226834/sampling-from-a-mixture-of-two-gamma-distributions/226837#226837
+ // this approach seems to be around 20000 times faster
+
+ // if one of the diffusions is off, the tail is not present
+ if (m_longitudinalDiffusionSigma == 0 || m_transverseDiffusionSigma == 0) {
+ float tmp = m_random->Gaus(0.0, posY*m_transverseDiffusionSigma);
+ // limit random number to be -1 < x < 1
+ while (std::abs(tmp) > 1.) tmp = m_random->Gaus(0.0, posY*m_transverseDiffusionSigma);
+ return tmp;
+ }
+ float tmp = generateTransverseDiffusion(posY);
+ while (std::abs(tmp) > 1.) {tmp = generateTransverseDiffusion(posY);}
+
+ return tmp;
}
-float MM_StripsResponseSimulation::getLongitudinalDiffusion(float posY){
- m_longitudinalDiffusionFunction->SetParameters(1.0, 0., posY*m_longitudinalDiffusionSigma);
- return m_longitudinalDiffusionFunction->GetRandom();
+float MM_StripsResponseSimulation::getLongitudinalDiffusion(float posY) {
+ float tmp = m_random->Gaus(0.0, posY*m_longitudinalDiffusionSigma);
+ // We only want random numbers between -5 and 5
+ while (std::abs(tmp) > 5) { tmp = m_random->Gaus(0.0, posY*m_longitudinalDiffusionSigma); }
+ return tmp;
}
-float MM_StripsResponseSimulation::getEffectiveCharge(){return m_polyaFunction->GetRandom();}
+float MM_StripsResponseSimulation::getEffectiveCharge() {return m_polyaFunction->GetRandom();}
+
+float MM_StripsResponseSimulation::getPathLengthTraveled() {
+ float rndGaus = m_random->Gaus(m_interactionDensityMean, m_interactionDensitySigma);
+
+
+ // gaussian random number should be in the range from 0 to 10
+ while (rndGaus < 0. || rndGaus > 10.) {rndGaus = m_random->Gaus(m_interactionDensityMean, m_interactionDensitySigma);}
-float MM_StripsResponseSimulation::getPathLengthTraveled(){
- return ( 1. / m_interactionDensityFunction->GetRandom() ) * -1. * log( m_random->Uniform() );
+ return ( 1. / rndGaus) * -1. * log( m_random->Uniform() );
}
--
GitLab