12 #include <klm/simulation/ScintillatorSimulator.h>
15 #include <klm/bklm/geometry/GeometryPar.h>
16 #include <klm/eklm/geometry/GeometryData.h>
19 #include <framework/dataobjects/EventMetaData.h>
20 #include <framework/datastore/StoreObjPtr.h>
21 #include <framework/gearbox/Unit.h>
22 #include <framework/logging/Logger.h>
35 static const char MemErr[] =
"Memory allocation error.";
74 double digitizationInitialTime,
bool debug) :
77 m_DigitizationInitialTime(digitizationInitialTime),
79 m_FPGAStat(c_ScintillatorFirmwareNoSignal),
87 double time, attenuationTime;
117 time = digPar->getADCSamplingTime() * i;
119 exp(-digPar->getPEAttenuationFrequency() * time) * attenuationTime /
120 digPar->getADCSamplingTime();
135 free(m_amplitudeDirect);
136 free(m_amplitudeReflected);
138 free(m_ADCAmplitude);
139 free(m_SignalTimeDependence);
140 free(m_SignalTimeDependenceDiff);
141 free(m_Photoelectrons);
142 free(m_PhotoelectronIndex);
143 free(m_PhotoelectronIndex2);
160 for (
int i = 0; i < m_DigPar->getNDigitizations(); i++) {
162 m_amplitudeDirect[i] = 0;
163 m_amplitudeReflected[i] = 0;
170 const std::multimap<uint16_t, const BKLMSimHit*>::iterator& firstHit,
171 const std::multimap<uint16_t, const BKLMSimHit*>::iterator& end)
173 m_stripName =
"strip_" + std::to_string(firstHit->first);
178 geoPar->
findModule(hit->getSection(), hit->getSector(), hit->getLayer());
180 2.0 * (hit->isPhiReadout() ?
181 module->getPhiScintHalfLength(hit->getStrip()) :
182 module->getZScintHalfLength(hit->getStrip()));
183 std::vector<const BKLMSimHit*> hits;
184 for (std::multimap<uint16_t, const BKLMSimHit*>::iterator it = firstHit;
186 hits.push_back(it->second);
187 std::sort(hits.begin(), hits.end(), compareBKLMSimHits);
188 for (std::vector<const BKLMSimHit*>::iterator it = hits.begin();
189 it != hits.end(); ++it) {
191 m_Energy = m_Energy + hit->getEnergyDeposit();
193 double nPhotons = hit->getEnergyDeposit() * m_DigPar->getNPEperMeV();
195 double sipmDistance = hit->getPropagationTime() *
196 m_DigPar->getFiberLightSpeed();
197 double time = hit->getTime() + hit->getPropagationTime();
199 m_MCTime = hit->getTime();
202 if (hit->getTime() < m_MCTime)
203 m_MCTime = hit->getTime();
204 if (time < m_SiPMMCTime)
207 int generatedPhotons = gRandom->Poisson(nPhotons);
208 generatePhotoelectrons(stripLength, sipmDistance, generatedPhotons,
209 hit->getTime(),
false);
210 if (m_DigPar->getMirrorReflectiveIndex() > 0) {
211 generatedPhotons = gRandom->Poisson(nPhotons);
212 generatePhotoelectrons(stripLength, sipmDistance, generatedPhotons,
213 hit->getTime(),
true);
220 const std::multimap<uint16_t, const EKLMSimHit*>::iterator& firstHit,
221 const std::multimap<uint16_t, const EKLMSimHit*>::iterator& end)
223 m_stripName =
"strip_" + std::to_string(firstHit->first);
227 hit->getStrip()) / CLHEP::mm *
Unit::mm;
228 std::vector<const EKLMSimHit*> hits;
229 for (std::multimap<uint16_t, const EKLMSimHit*>::iterator it = firstHit;
231 hits.push_back(it->second);
232 std::sort(hits.begin(), hits.end(), compareEKLMSimHits);
233 for (std::vector<const EKLMSimHit*>::iterator it = hits.begin();
234 it != hits.end(); ++it) {
236 m_Energy = m_Energy + hit->getEnergyDeposit();
238 double nPhotons = hit->getEnergyDeposit() * m_DigPar->getNPEperMeV();
240 double sipmDistance = 0.5 * stripLength - hit->getLocalPosition().x();
241 double time = hit->getTime() +
242 sipmDistance / m_DigPar->getFiberLightSpeed();
246 m_MCTime = time < m_MCTime ? time : m_MCTime;
247 int generatedPhotons = gRandom->Poisson(nPhotons);
248 generatePhotoelectrons(stripLength, sipmDistance, generatedPhotons,
249 hit->getTime(),
false);
250 if (m_DigPar->getMirrorReflectiveIndex() > 0) {
251 generatedPhotons = gRandom->Poisson(nPhotons);
252 generatePhotoelectrons(stripLength, sipmDistance, generatedPhotons,
253 hit->getTime(),
true);
262 fillSiPMOutput(m_amplitudeDirect,
true,
false);
263 fillSiPMOutput(m_amplitudeReflected,
false,
true);
264 for (
int i = 0; i < m_DigPar->getNDigitizations(); i++)
265 m_amplitude[i] = m_amplitudeDirect[i] + m_amplitudeReflected[i];
267 fillSiPMOutput(m_amplitude,
true,
true);
269 if (m_DigPar->getMeanSiPMNoise() > 0)
270 addRandomSiPMNoise();
272 m_FPGAStat = m_fitter->fit(m_ADCAmplitude, m_Threshold, &m_FPGAFit);
273 if (m_FPGAStat != c_ScintillatorFirmwareSuccessfulFit)
284 for (i = 0; i < m_DigPar->getNDigitizations(); i++)
285 m_amplitude[i] = m_amplitude[i] +
286 gRandom->Poisson(m_DigPar->getMeanSiPMNoise());
291 int* currentIndexArray, *newIndexArray, *tmpIndexArray;
292 int i, i1, i2, i1Max, i2Max, j, mergeSize;
293 currentIndexArray = m_PhotoelectronIndex;
294 newIndexArray = m_PhotoelectronIndex2;
296 while (mergeSize < nPhotoelectrons) {
297 for (i = 0; i < nPhotoelectrons; i = i + 2 * mergeSize) {
301 if (i2 > nPhotoelectrons)
302 i2 = nPhotoelectrons;
304 i2Max = i2 + mergeSize;
305 if (i2Max > nPhotoelectrons)
306 i2Max = nPhotoelectrons;
307 while (i1 < i1Max || i2 < i2Max) {
310 if (m_Photoelectrons[currentIndexArray[i1]].bin <
311 m_Photoelectrons[currentIndexArray[i2]].bin) {
312 newIndexArray[j] = currentIndexArray[i1];
315 newIndexArray[j] = currentIndexArray[i2];
319 newIndexArray[j] = currentIndexArray[i1];
323 newIndexArray[j] = currentIndexArray[i2];
329 tmpIndexArray = currentIndexArray;
330 currentIndexArray = newIndexArray;
331 newIndexArray = tmpIndexArray;
332 mergeSize = mergeSize * 2;
334 return currentIndexArray;
338 double stripLen,
double distSiPM,
int nPhotons,
double timeShift,
341 const double maxHitTime = m_DigPar->getNDigitizations() *
342 m_DigPar->getADCSamplingTime();
345 double hitTime, deExcitationTime, cosTheta, hitDist, selection;
346 double inverseLightSpeed, inverseAttenuationLength;
347 inverseLightSpeed = 1.0 / m_DigPar->getFiberLightSpeed();
348 inverseAttenuationLength = 1.0 / m_DigPar->getAttenuationLength();
350 for (i = 0; i < nPhotons; i++) {
351 if (m_npe >= m_PhotoelectronBufferSize)
352 reallocPhotoElectronBuffers(m_npe + 100);
353 cosTheta = gRandom->Uniform(m_DigPar->getMinCosTheta(), 1);
355 hitDist = distSiPM / cosTheta;
357 hitDist = (2.0 * stripLen - distSiPM) / cosTheta;
359 selection = gRandom->Uniform();
360 if (selection > exp(-hitDist * inverseAttenuationLength))
364 selection = gRandom->Uniform();
365 if (selection > m_DigPar->getMirrorReflectiveIndex())
369 gRandom->Exp(m_DigPar->getScintillatorDeExcitationTime()) +
370 gRandom->Exp(m_DigPar->getFiberDeExcitationTime());
371 hitTime = hitDist * inverseLightSpeed + deExcitationTime +
372 timeShift - m_DigitizationInitialTime;
373 if (hitTime >= maxHitTime)
376 m_Photoelectrons[m_npe].bin =
377 floor(hitTime / m_DigPar->getADCSamplingTime());
379 m_Photoelectrons[m_npe].bin = -1;
380 m_Photoelectrons[m_npe].expTime =
381 exp(m_DigPar->getPEAttenuationFrequency() * hitTime);
382 m_Photoelectrons[m_npe].isReflected = isReflected;
383 m_PhotoelectronIndex[m_npe] = m_npe;
451 double attenuationTime, sig, expSum;
453 int ind1, ind2, ind3;
457 attenuationTime = 1.0 / m_DigPar->getPEAttenuationFrequency();
458 indexArray = sortPhotoelectrons(m_npe);
463 bin = m_Photoelectrons[indexArray[ind1]].bin;
468 if (bin != m_Photoelectrons[indexArray[ind2]].bin)
472 for (ind3 = ind1; ind3 != ind2; ind3++) {
473 if (m_Photoelectrons[indexArray[ind3]].isReflected && !useReflected)
475 if (!m_Photoelectrons[indexArray[ind3]].isReflected && !useDirect)
478 sig = attenuationTime - m_Photoelectrons[indexArray[ind3]].expTime *
479 m_SignalTimeDependence[bin + 1];
480 hist[bin] = hist[bin] + sig;
482 expSum = expSum + m_Photoelectrons[indexArray[ind3]].expTime;
485 maxBin = m_DigPar->getNDigitizations() - 1;
487 maxBin = m_Photoelectrons[indexArray[ind2]].bin;
488 for (i = bin + 1; i <= maxBin; i++) {
489 sig = m_SignalTimeDependenceDiff[i] * expSum;
490 hist[i] = hist[i] + sig;
503 if (m_Pedestal == 0 || m_PhotoelectronAmplitude == 0)
504 B2FATAL(
"Incorrect EKLM ADC simulation parameters.");
505 for (i = 0; i < m_DigPar->getNDigitizations(); i++) {
506 amp = m_Pedestal - m_PhotoelectronAmplitude * m_amplitude[i];
507 if (amp < m_DigPar->getADCSaturation())
508 amp = m_DigPar->getADCSaturation();
509 m_ADCAmplitude[i] = floor(amp);
526 intg = m_FPGAFit.getAmplitude();
527 return intg * m_DigPar->getPEAttenuationFrequency() /
528 m_PhotoelectronAmplitude;
546 TFile* hfile =
nullptr;
547 TH1D* histAmplitudeDirect =
nullptr;
548 TH1D* histAmplitudeReflected =
nullptr;
549 TH1D* histAmplitude =
nullptr;
550 TH1D* histADCAmplitude =
nullptr;
552 histAmplitudeDirect =
553 new TH1D(
"histAmplitudeDirect", m_stripName.c_str(),
554 m_DigPar->getNDigitizations(), 0, m_histRange);
555 histAmplitudeReflected =
556 new TH1D(
"histAmplitudeReflected", m_stripName.c_str(),
557 m_DigPar->getNDigitizations(), 0, m_histRange);
559 new TH1D(
"histAmplitude", m_stripName.c_str(),
560 m_DigPar->getNDigitizations(), 0, m_histRange);
562 new TH1D(
"histADCAmplitude", m_stripName.c_str(),
563 m_DigPar->getNDigitizations(), 0, m_histRange);
564 }
catch (std::bad_alloc& ba) {
567 for (i = 0; i < m_DigPar->getNDigitizations(); i++) {
568 histAmplitudeDirect->SetBinContent(i + 1, m_amplitudeDirect[i]);
569 histAmplitudeReflected->SetBinContent(i + 1, m_amplitudeReflected[i]);
570 histAmplitude->SetBinContent(i + 1, m_amplitude[i]);
571 histADCAmplitude->SetBinContent(i + 1, m_ADCAmplitude[i]);
573 str = std::string(
"experiment_") + std::to_string(event->getExperiment()) +
574 "_run_" + std::to_string(event->getRun()) +
"_event_" +
575 std::to_string(event->getEvent()) +
"_" + m_stripName +
".root";
577 hfile =
new TFile(str.c_str(),
"NEW");
578 }
catch (std::bad_alloc& ba) {
581 hfile->Append(histAmplitudeDirect);
582 hfile->Append(histAmplitudeReflected);
583 hfile->Append(histAmplitude);
584 hfile->Append(histADCAmplitude);