Merge branch 'vmmSpeedupRound2_master' into 'master'

Manual sweep of "Speedup of the MM electronics simulation "

See merge request atlas/athena!38020

MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_ElectronicsToolInput.h

@@ -24,9 +24,9 @@ class MM_ElectronicsToolInput {
 
     ~MM_ElectronicsToolInput() {}
 
-    std::vector<int> NumberOfStripsPos() const { return m_NumberOfStripsPos; }
-    std::vector<std::vector<float>> chipCharge() const { return m_chipCharge; }
-    std::vector<std::vector<float>> chipTime() const { return m_chipTime; }
+   const std::vector<int>& NumberOfStripsPos() const { return m_NumberOfStripsPos; }
+   const std::vector<std::vector<float>>& chipCharge() const { return m_chipCharge; }
+   const std::vector<std::vector<float>>& chipTime() const { return m_chipTime; }
     Identifier digitID() const { return m_digitID; }
     float kineticEnergy() const { return m_kineticEnergy; }
 

MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripsResponseSimulation.h

@@ -181,5 +181,13 @@ private:
   //Declaring private message stream member.
   mutable Athena::MsgStreamMember m_msg = Athena::MsgStreamMember("MMStripResponseSimulation");
 
+  // seperate random number generation for performance monitoring
+  float getTransversDiffusion(float posY);
+  float getLongitudinalDiffusion(float posY);
+  float getEffectiveCharge();
+  float getPathLengthTraveled();
+
+
+
 };
 #endif

MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/VMM_Shaper.h

@@ -37,8 +37,9 @@ class VMM_Shaper{
      double vmmResponse(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, double time) const;
      double findPeak(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, const double electronicsThreshold) const;
      bool aboveThresholdSimple(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, const double electronicsThreshold) const;
-
-
+     double m_pole0_ns;
+     double m_re_pole1_ns; 
+     double m_im_pole1_ns;
 };
 
 #endif  // MM_DIGITIZATION_VMM_SHAPER_H

MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripsResponseSimulation.cxx

@@ -114,7 +114,11 @@ void MM_StripsResponseSimulation::initFunctions()
 
 	m_longitudinalDiffusionFunction = new TF1("longdiff","gaus", -5., 5.);
 
-	m_transverseDiffusionFunction = new TF1("transdiff", "1.*TMath::Exp(-TMath::Power(x,2.)/(2.*[0]*[0])) + 0.001*TMath::Exp(-TMath::Power(x,2)/(2.*[1]*[1]))", -1., 1.);
+	//m_transverseDiffusionFunction = new TF1("transdiff", "1.*TMath::Exp(-TMath::Power(x,2.)/(2.*[0]*[0])) + 0.001*TMath::Exp(-TMath::Power(x,2)/(2.*[1]*[1]))", -1., 1.);
+	m_transverseDiffusionFunction = new TF1("transdiff", "gaus(0) + gaus(3)", -1., 1.);
+  m_transverseDiffusionFunction->SetParameters(1.,0.,1.,0.001,0.,1.);
+	//m_transverseDiffusionFunction = new TF1("transdiff", "gaus", -1., 1.);
+  //m_transverseDiffusionFunction->SetParameters(1.,0.,1.);
 
 	m_random = new TRandom3(0);
 
@@ -202,7 +206,7 @@ void MM_StripsResponseSimulation::whichStrips( const float & hitx,
 	ATH_MSG_DEBUG("LorentzAngle vs theta: " <<lorentzAngle <<" " <<theta);
     ATH_MSG_DEBUG("Function pointer points to " << m_interactionDensityFunction);
 
-	float pathLengthTraveled = ( 1. / m_interactionDensityFunction->GetRandom() ) * -1. * std::log( m_random->Uniform() );
+	float pathLengthTraveled = getPathLengthTraveled(); 
 
 	float pathLength  = m_driftGapWidth/std::cos(theta); // Increasing path length based on XZ angle
 	pathLength       /= std::cos(alpha);                 // Further increasing path length for YZ angle
@@ -231,19 +235,19 @@ void MM_StripsResponseSimulation::whichStrips( const float & hitx,
 
 		for (auto& Electron : IonizationCluster.getElectrons()){
 
-			m_longitudinalDiffusionFunction->SetParameters(1.0, 0., initialPosition.Y()*m_longitudinalDiffusionSigma);
+		//	m_longitudinalDiffusionFunction->SetParameters(1.0, 0., initialPosition.Y()*m_longitudinalDiffusionSigma);
 
-			if ( m_longitudinalDiffusionSigma == 0 || m_transverseDiffusionSigma == 0) {
+		//	if ( m_longitudinalDiffusionSigma == 0 || m_transverseDiffusionSigma == 0) {
 
-				m_transverseDiffusionFunction->SetParameters( initialPosition.Y()*m_transverseDiffusionSigma , 0.0);
+		//		m_transverseDiffusionFunction->SetParameters( initialPosition.Y()*m_transverseDiffusionSigma , 0.0);
 
-			} else {
+		//	} else {
 
-				m_transverseDiffusionFunction->SetParameters( initialPosition.Y()*m_transverseDiffusionSigma , 1.0);
+		//		m_transverseDiffusionFunction->SetParameters( initialPosition.Y()*m_transverseDiffusionSigma , 1.0);
 
-			}
+		//	}
 
-			Electron->setOffsetPosition(m_transverseDiffusionFunction->GetRandom(), m_longitudinalDiffusionFunction->GetRandom());
+			Electron->setOffsetPosition(getTransversDiffusion(initialPosition.Y()) , getLongitudinalDiffusion(initialPosition.Y()) );
 
 		}
 
@@ -254,7 +258,8 @@ void MM_StripsResponseSimulation::whichStrips( const float & hitx,
 
 		for (auto& Electron : IonizationCluster.getElectrons()){
 
-			float effectiveCharge = m_polyaFunction->GetRandom();
+			//float effectiveCharge = m_polyaFunction->GetRandom();
+			float effectiveCharge = getEffectiveCharge();
 
 			Electron->setCharge( effectiveCharge );
 
@@ -275,7 +280,7 @@ void MM_StripsResponseSimulation::whichStrips( const float & hitx,
 		//---
 		m_IonizationClusters.push_back(IonizationCluster);
 
-		pathLengthTraveled +=  (1. / m_interactionDensityFunction->GetRandom() ) * -1. * std::log( m_random->Uniform() );
+		pathLengthTraveled +=  getPathLengthTraveled();
 
 		ATH_MSG_DEBUG("Path length traveled: " << pathLengthTraveled);
 
@@ -332,6 +337,33 @@ 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::getLongitudinalDiffusion(float posY){
+  m_longitudinalDiffusionFunction->SetParameters(1.0, 0., posY*m_longitudinalDiffusionSigma);
+  return m_longitudinalDiffusionFunction->GetRandom();
+}
+
+float  MM_StripsResponseSimulation::getEffectiveCharge(){return m_polyaFunction->GetRandom();}
+
+float MM_StripsResponseSimulation::getPathLengthTraveled(){
+	return  ( 1. / m_interactionDensityFunction->GetRandom() ) * -1. * log( m_random->Uniform() );
+}
+
+
+
+
 /*******************************************************************************/
 MM_StripsResponseSimulation::~MM_StripsResponseSimulation()
 {

MuonSpectrometer/MuonDigitization/MM_Digitization/src/VMM_Shaper.cxx

@@ -42,16 +42,20 @@ void VMM_Shaper::initialize() {
     
     // preCalculate factor to avoid recalculating for each electron
     m_preCalculationVMMShaper = chargeScaleFactor*m_peakTimeChargeScaling*std::pow(m_a, 3)*m_pole0*m_pole1_square;
+
+    m_pole0_ns = m_pole0*(10^-9);
+    m_re_pole1_ns = m_re_pole1*(10^-9);
+    m_im_pole1_ns = m_im_pole1*(10^-9);
 }
 
 double VMM_Shaper::vmmResponse(const std::vector<float> &effectiveCharge, const std::vector<float> &electronsTime, double time) const{
     double response = 0;
     for (unsigned int i_electron = 0; i_electron < effectiveCharge.size(); i_electron++) {
             if (time < electronsTime.at(i_electron)) continue;
-            double t = (time-electronsTime.at(i_electron))*(10^-9);
+            double t = (time-electronsTime.at(i_electron));
             // now follows the vmm shaper response function provided by G. Iakovidis
             // It is described in section 7.1.3 of https://cds.cern.ch/record/1955475
-            double st = effectiveCharge.at(i_electron)*m_preCalculationVMMShaper*((K0*std::exp(-t*m_pole0))+(2.*m_k1_abs*std::exp(-t*m_re_pole1)*std::cos(-t*m_im_pole1+m_argK1)));
+            double st = effectiveCharge.at(i_electron)*m_preCalculationVMMShaper*((K0*std::exp(-t*m_pole0_ns))+(2.*m_k1_abs*std::exp(-t*m_re_pole1_ns)*std::cos(-t*m_im_pole1_ns+m_argK1)));
             response += st;
     }
     return response;
@@ -74,13 +78,29 @@ void VMM_Shaper::vmmThresholdResponse(const std::vector<float> &effectiveCharge,
     double startTime = m_lowerTimeWindow;
     double minElectronTime = *std::min_element(electronsTime.begin(), electronsTime.end());
     if (startTime < minElectronTime) startTime = minElectronTime;  // if smallest strip times are higher then the lower time window, just start the loop from the smallest electron time
-
+    
     double tmpTimeAtThreshold = -9999;
-    for (double time = startTime; time < m_upperTimeWindow; time += m_timeStep) {
-        if (vmmResponse(effectiveCharge, electronsTime, time) >= electronicsThreshold) {
-            tmpTimeAtThreshold = time;
-            break;
-        }
+    
+    for (double time = startTime; time < minElectronTime + 0.9*m_peakTime; time += 10*m_timeStep) { // quick search till the first possible peak
+      if(vmmResponse(effectiveCharge, electronsTime, time) <= electronicsThreshold ) continue;
+      for(double fineTime = time; fineTime >= time - 10*m_timeStep; fineTime -= m_timeStep) { // since value above threshold was found, loop back in time to find the crossing with the timeStep precission
+          if(vmmResponse(effectiveCharge, electronsTime, fineTime) >= electronicsThreshold) continue;
+              tmpTimeAtThreshold = fineTime + 0.5*m_timeStep; //  get time between time above and time below threshold 
+              break;
+          }
+      break;
+    } 
+    
+    if(tmpTimeAtThreshold == -9999) {  // threshold crossing not yet found 
+        // check if first possible peak was before start time
+        double tmpStartTime = std::max(minElectronTime + 0.9*m_peakTime, startTime); 
+        for (double time = tmpStartTime; time < m_upperTimeWindow; time += m_timeStep) {
+          if (vmmResponse(effectiveCharge, electronsTime, time) >= electronicsThreshold) {
+              tmpTimeAtThreshold = time;
+              break;
+          }
+       } 
+    
     }
 
     if (tmpTimeAtThreshold == -9999) return;
@@ -97,44 +117,58 @@ double VMM_Shaper::findPeak(const std::vector<float> &effectiveCharge, const std
     if(effectiveCharge.size()==0) return -9999; // protect min_element
     double startTime = m_lowerTimeWindow;
     double minElectronTime = *std::min_element(electronsTime.begin(), electronsTime.end());
+
+    minElectronTime += 0.8*m_peakTime;  // only start looking for the peak close to the first possible peak
     if(startTime < minElectronTime) startTime = minElectronTime; // if smallest strip times are higher then the lower time window, just start the loop from the smallest electron time
 
     double oldResponse = 0;
-    double oldDerivative = 0 , currentDerivative = 0;
+    //double currentDerivative = 0;
 
-    for (double time = startTime; time < m_upperTimeWindow; time += m_timeStep) {
+    double timeStepScaleFactor = 5.0;
 
-        double response = vmmResponse(effectiveCharge, electronsTime, time);
+    for (double time = startTime; time < m_upperTimeWindow; time += m_timeStep*timeStepScaleFactor) {
 
-        oldDerivative = currentDerivative;
-        currentDerivative = (response-oldResponse) * m_inverseTimeStep;
-        //  check if sign of derivative has not changed ==> no peak;  or if response is below threshold
-        if (oldDerivative*currentDerivative >= 0 || oldResponse < electronicsThreshold) {oldResponse = response; continue;}
+        double response = vmmResponse(effectiveCharge, electronsTime, time);
+        if(oldResponse < response){oldResponse=response; continue;}
+        oldResponse = response;
 
-        // from here one its assumed that a peak above threshold was found
+        int searchWindow = 5;
 
-        //  check if the derivative is monoton falling within the given window of 5 bins
-        bool checkDerivative = true;
-        double tmp_checkOldDerivative = 0, tmp_checkCurrentDerivative = 0;
-        double tmp_checkOldResponse = 0;
+        std::vector<double> tmpTime, tmpResponse;
 
-        int searchWindow = 5;
+        tmpTime.reserve(2*timeStepScaleFactor);
+        tmpResponse.reserve(2*timeStepScaleFactor);
 
-        for (int i_timeOfPeak = -searchWindow; i_timeOfPeak <= searchWindow; i_timeOfPeak ++) {
-           double response = vmmResponse(effectiveCharge, electronsTime, time + i_timeOfPeak * m_timeStep);
 
-           tmp_checkOldDerivative = tmp_checkCurrentDerivative;
-           tmp_checkCurrentDerivative = (response - tmp_checkOldResponse) * m_inverseTimeStep;
-           tmp_checkOldResponse = response;
+        for(double fineTime = (time-1.5*m_timeStep*timeStepScaleFactor); fineTime < time+0.5*m_timeStep*timeStepScaleFactor; fineTime += m_timeStep) {
+          tmpTime.push_back(fineTime);
+          tmpResponse.push_back(vmmResponse(effectiveCharge,electronsTime,fineTime));     
+        }
 
-           if (i_timeOfPeak >= -searchWindow + 2  // needs two iterations to fill the variables
-                && tmp_checkOldDerivative < tmp_checkCurrentDerivative) {  // derivative is not falling monothonic
-               checkDerivative = false;
-               break;
+        int nBins = tmpTime.size();
+
+        for(int i_time = 1; i_time < nBins-1;i_time++){
+            if(tmpResponse.at(i_time)<tmpResponse.at(i_time+1)) continue;
+
+            if(tmpResponse.at(i_time) < electronicsThreshold) break; 
+            
+            bool checkTimeWindow = false;
+            for(int i_timeOfPeak = i_time - searchWindow + 1; i_timeOfPeak <= i_time + searchWindow; i_timeOfPeak++) {
+              if(i_timeOfPeak < 1 || i_timeOfPeak == nBins - 1) continue;
+              
+              double oldDerivative = (tmpResponse.at(i_time) - tmpResponse.at(i_time-1));
+              double newDerivative = (tmpResponse.at(i_time+1) - tmpResponse.at(i_time));
+              if (newDerivative > oldDerivative){
+                checkTimeWindow = false;
+                break; 
+              } else {
+                checkTimeWindow = true;
+              }
+            
            }
+           if(checkTimeWindow) return tmpTime.at(i_time);
+
         }
-        if (!checkDerivative) continue;
-        return time - m_timeStep;
     }
     return -9999;  // no peak found
 }