diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_ElectronicsToolInput.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_ElectronicsToolInput.h
index b818d16c7ea684ce7b4e1d8ee4e2d1dd54ce27a3..62683cca8507063d668ffaac19f30ae3fd56d673 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_ElectronicsToolInput.h
+++ b/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; }
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripsResponseSimulation.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripsResponseSimulation.h
index b879888dabe75e8682b830cdee757ff36533ab43..2ed48255b45b5c69f12b6fed4a51d221e8c11495 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/MM_StripsResponseSimulation.h
+++ b/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
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/VMM_Shaper.h b/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/VMM_Shaper.h
index 4def653f55ca39a53c6d4ff82f3b8fe8a1e38a9f..8b4a6672c970cd217a60e7c777af46e2e8acea9b 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/MM_Digitization/VMM_Shaper.h
+++ b/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
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripsResponseSimulation.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripsResponseSimulation.cxx
index 2e71f603fb9fc3bc183e9f498273aed4d4835e63..18216e7d3c9db0b79416d9b7b12e8544da4b8b83 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/MM_StripsResponseSimulation.cxx
+++ b/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()
{
diff --git a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/VMM_Shaper.cxx b/MuonSpectrometer/MuonDigitization/MM_Digitization/src/VMM_Shaper.cxx
index 39d0763c2ab71b15c51ad1432226b15fc35274cc..788f87ecec6b6e75b8c89ea2cac3c4f98dcf37a4 100644
--- a/MuonSpectrometer/MuonDigitization/MM_Digitization/src/VMM_Shaper.cxx
+++ b/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
}