15 #include "CLHEP/Random/RandFlat.h"
16 #include "CLHEP/Random/RandGamma.h"
17 #include "CLHEP/Random/RandGaussZiggurat.h"
18 #include "CLHEP/Random/RandGeneral.h"
19 #include "CLHEP/Random/RandomEngine.h"
20 #include "GaudiKernel/SystemOfUnits.h"
34 constexpr
float electronsToFC = 1. / 6241.;
35 constexpr
float polyaTheta = 1.765;
42 m_cfg{std::move(
mod)} {
47 std::vector<double> NelectronProb{65.6, 15., 6.4, 3.5, 2.25, 1.55, 1.05, 0.81, 0.61, 0.49, 0.39, 0.3,
48 0.25, 0.2, 0.16, 0.12, 0.095, 0.075, 0.063, 0.054, 0.049, 0.045, 0.044};
49 NelectronProb.reserve(s_NelectronPropBins);
50 for (
unsigned int Nelectron = NelectronProb.size(); Nelectron < s_NelectronPropBins; ++Nelectron)
51 NelectronProb.push_back(21.6 / ((Nelectron) * (Nelectron)));
53 m_randNelectrons = std::make_unique<CLHEP::RandGeneral>(NelectronProb.data(), NelectronProb.size(), 1);
55 ATH_MSG_DEBUG(
"MM_StripsResponseSimulation::initializationFrom set values");
61 m_outputFile = std::make_unique<TFile>(
"MM_StripsResponse_Plots.root",
"RECREATE");
64 m_mapOfHistograms[
"nInteractions"] = std::make_unique<TH1F>(
"nInteractions",
"nInteractions", 200, 0, 200);
65 m_mapOfHistograms[
"nElectrons"] = std::make_unique<TH1F>(
"nElectrons",
"nElectrons", 200, 0, 200);
66 m_mapOfHistograms[
"lorentzAngle"] = std::make_unique<TH1F>(
"lorentzAngle",
"lorentzAngle", 100, 0, 3);
67 m_mapOfHistograms[
"effectiveCharge"] = std::make_unique<TH1F>(
"effectiveCharge",
"effectiveCharge", 100, 0, 1e5);
68 m_mapOfHistograms[
"effectiveNElectrons"] = std::make_unique<TH1F>(
"effectiveNElectrons",
"effectiveNElectrons", 200, 0, 200);
69 m_mapOfHistograms[
"totalEffectiveCharge"] = std::make_unique<TH1F>(
"totalEffectiveCharge",
"totalEffectiveCharge", 400, 0, 200);
70 m_mapOf2DHistograms[
"totalEffectiveChargeVsTheta"] = std::make_unique<TH2F>(
"totalEffectiveChargeVsTheta",
"totalEffectiveChargeVsTheta", 17, 0, 35, 300, 0, 300);
72 m_mapOf2DHistograms[
"pathLengthVsTheta"] = std::make_unique<TH2F>(
"pathLengthVsTheta",
"pathLengthVsTheta", 100, 0, 10, 100, -3, 3);
73 m_mapOf2DHistograms[
"pathLengthVsAlpha"] = std::make_unique<TH2F>(
"pathLengthVsAlpha",
"pathLengthVsAlpha", 100, 0, 10, 100, -3, 3);
75 m_mapOf2DHistograms[
"lorentzAngleVsTheta"] = std::make_unique<TH2F>(
"lorentzAngleVsTheta",
"lorentzAngleVsTheta", 100, -3, 3, 100, -3, 3);
76 m_mapOf2DHistograms[
"lorentzAngleVsBy"] = std::make_unique<TH2F>(
"lorentzAngleVsBy",
"lorentzAngleVsBy", 100, -3, 3, 100, -2., 2.);
79 std::make_unique<TH2F>(
"deltaLorentzAngleVsByForPosTheta",
"deltaLorentzAngleVsByForPosTheta", 100, -3, 3, 100, -2., 2.);
81 std::make_unique<TH2F>(
"deltaLorentzAngleVsByForNegTheta",
"deltaLorentzAngleVsByForNegTheta", 100, -3, 3, 100, -2., 2.);
84 std::make_unique<TH2F>(
"driftDistanceVsDriftTime",
"driftDistanceVsDriftTime", 100, 0, 10, 100, 0, 200);
103 double gainFraction,
double stripPitch,
104 CLHEP::HepRandomEngine* rndmEngine)
const {
107 whichStrips(cache, digiInput, gainFraction, stripPitch, rndmEngine);
112 std::move(cache.finalqStrip),
113 std::move(cache.finaltStrip));
115 return tmpStripToolOutput;
124 double gainFraction,
double stripPitch,
125 CLHEP::HepRandomEngine* rndmEngine)
const {
126 ATH_MSG_DEBUG(
"Starting to calculate strips that got fired");
131 const int stripMinID = digiInput.
stripMinID();
132 const int stripMaxID = digiInput.
stripMaxID();
135 const float eventTime = digiInput.
eventTime();
139 int nPrimaryIons = 0;
147 std::lock_guard guard{fillMutex};
160 std::lock_guard guard{fillMutex};
169 while (pathLengthTraveled < pathLength) {
171 std::unique_ptr<MM_IonizationCluster> IonizationCluster =
180 ATH_MSG_DEBUG(
"New interaction starting at x,y, pathLengthTraveled: " << initialPosition.x() <<
" " << initialPosition.y() <<
" "
181 << pathLengthTraveled);
190 int tmpEffectiveNElectrons = 0;
195 Electron->setCharge(effectiveCharge);
197 ATH_MSG_DEBUG(
"Electron's effective charge is " << effectiveCharge);
200 std::lock_guard guard{fillMutex};
207 tmpEffectiveNElectrons += effectiveCharge;
211 std::lock_guard guard{fillMutex};
219 ATH_MSG_DEBUG(
"Path length traveled: " << pathLengthTraveled);
226 float timeresolution = 0.01;
245 std::lock_guard guard{fillMutex};
250 stripResponseObject.
totalCharge() * electronsToFC);
256 std::lock_guard guard{fillMutex};
257 static std::atomic<unsigned> written{0};
259 for (
size_t iIonization = 0; iIonization < cache.
IonizationClusters.size(); ++iIonization) {
261 grIonizationXZ.SetPoint(iIonization, ionizationPosition.x(), ionizationPosition.y());
263 m_outputFile->WriteObject(&grIonizationXZ, Form(
"ionizationXZ_%u", ++written));
271 m_outputFile->WriteObject(&grElectronsXZ, Form(
"electronsXZ_%u", ++written));
286 const double uni = CLHEP::RandFlat::shoot(rndmEngine, 0, 1.0 +
scale);
287 if (uni <
scale)
return CLHEP::RandGaussZiggurat::shoot(rndmEngine, 0.0, 1.0);
307 }
while (std::abs(diffusion) > 1.);
312 ATH_MSG_VERBOSE(
"getTransverseDiffusion() -- posY: "<<
posY<<
" engine seed: "<<rndmEngine->getSeed()<<
" "<<rndmEngine->flat()<<
" diffusion "<<diffusion);
313 }
while (std::abs(diffusion)> 1.);
322 ATH_MSG_VERBOSE(
"getLongitudinalDiffusion() -- posY: "<<
posY<<
" engine seed: "<<rndmEngine->getSeed()<<
" "<<rndmEngine->flat()<<
" diffusion "<<diffusion);
323 }
while (std::abs(diffusion) > 5);
331 return CLHEP::RandGamma::shoot(rndmEngine, 1. + polyaTheta, (1. + polyaTheta) /
m_cfg.
avalancheGain);
341 ATH_MSG_VERBOSE(
"getPathLengthTraveled() -- engine seed: "<<rndmEngine->getSeed()<<
" "<<rndmEngine->flat()<<
" rndGaus: "<<rndGaus);
342 }
while (rndGaus < 0. || rndGaus > 10.) ;
344 return (1. / rndGaus) * -1. *
std::log(CLHEP::RandFlat::shoot(rndmEngine));