26#include <fairlogger/Logger.h>
38 LOG(info) <<
"Initializing digitizer";
40 mChips.resize(mNumberOfChips);
41 for (
int i = mNumberOfChips;
i--;) {
42 mChips[
i].setChipIndex(
i);
44 mChips[
i].setNoiseMap(mNoiseMap);
48 mChips[
i].setDeadChanMap(mDeadChanMap);
53 auto file = TFile::Open(mResponseFile.data());
55 LOG(fatal) <<
"Cannot open response file " << mResponseFile;
60 mChipSimRespVD = mChipSimResp;
61 mChipSimRespMLOT = mChipSimResp;
65 LOG(info) <<
" Depth max VD: " << mChipSimRespVD->getDepthMax();
66 LOG(info) <<
" Depth min VD: " << mChipSimRespVD->getDepthMin();
68 LOG(info) <<
" Depth max MLOT: " << mChipSimRespMLOT->getDepthMax();
69 LOG(info) <<
" Depth min MLOT: " << mChipSimRespMLOT->getDepthMin();
71 float thicknessVD = 0.0095;
72 float thicknessMLOT = 0.1;
76 mSimRespVDShift = -mChipSimRespVD->getDepthMax();
79 mSimRespMLOTShift = mChipSimRespMLOT->getDepthMax() - thicknessMLOT / 2.f;
80 mSimRespOrientation =
false;
85 LOGP(info,
"TRK Digitizer is initialised.");
87 LOGP(info,
"VD shift = {} ; ML/OT shift = {} = {} - {}", mSimRespVDShift, mSimRespMLOTShift, mChipSimRespMLOT->getDepthMax(), thicknessMLOT / 2.f);
88 LOGP(info,
"VD pixel scale on x = {} ; z = {}", mSimRespVDScaleX, mSimRespVDScaleZ);
89 LOGP(info,
"ML/OT pixel scale on x = {} ; z = {}", mSimRespMLOTScaleX, mSimRespMLOTScaleZ);
90 LOGP(info,
"Response orientation: {}", mSimRespOrientation ?
"flipped" :
"normal");
97 if (mGeometry->getSubDetID(chipID) == 0) {
98 return mChipSimRespVD;
101 else if (mGeometry->getSubDetID(chipID) == 1) {
102 return mChipSimRespMLOT;
112 LOG(info) <<
" Digitizing " << mGeometry->
getName() <<
" (ID: " << mGeometry->
getDetID()
113 <<
") hits of entry " << evID <<
" from source " << srcID
114 <<
" at time " << mEventTime <<
" ROFrame= " << mNewROFrame <<
")"
116 <<
" Min/Max ROFrames " << mROFrameMin <<
"/" << mROFrameMax;
122 if (mNewROFrame > mROFrameMin) {
126 int nHits = hits->size();
127 std::vector<int> hitIdx(nHits);
128 std::iota(std::begin(hitIdx), std::end(hitIdx), 0);
130 std::sort(hitIdx.begin(), hitIdx.end(),
131 [hits](
auto lhs,
auto rhs) {
132 return (*hits)[lhs].GetDetectorID() < (*hits)[rhs].GetDetectorID();
134 LOG(info) <<
"Processing " << nHits <<
" hits";
135 for (
int i : hitIdx) {
136 processHit((*hits)[
i], mROFrameMax, evID, srcID);
150 LOG(info) <<
"Setting event time ";
162 if (mCollisionTimeWrtROF < 0 && nbc > 0) {
168 LOG(info) <<
" NewROFrame " << mNewROFrame <<
" = " << nbc <<
"/" << mParams.
getROFrameLengthInBC() <<
" (nbc/mParams.getROFrameLengthInBC()";
176 if (mNewROFrame < mROFrameMin) {
177 LOG(error) <<
"New ROFrame " << mNewROFrame <<
" (" << irt <<
") precedes currently cashed " << mROFrameMin;
178 throw std::runtime_error(
"deduced ROFrame precedes already processed one");
181 if (mParams.
isContinuous() && mROFrameMax < mNewROFrame) {
182 mROFrameMax = mNewROFrame - 1;
190 if (frameLast > mROFrameMax) {
191 frameLast = mROFrameMax;
194 getExtraDigBuffer(mROFrameMax);
195 LOG(info) <<
"Filling " << mGeometry->
getName() <<
" digits output for RO frames " << mROFrameMin <<
":"
201 for (; mROFrameMin <= frameLast; mROFrameMin++) {
205 auto& extra = *(mExtraBuff.front().get());
206 for (
auto& chip : mChips) {
207 if (chip.isDisabled()) {
211 auto&
buffer = chip.getPreDigits();
215 auto itBeg =
buffer.begin();
217 ULong64_t maxKey = chip.getOrderingKey(mROFrameMin + 1, 0, 0) - 1;
218 for (; iter !=
buffer.end(); ++iter) {
219 if (iter->first > maxKey) {
222 auto& preDig = iter->second;
224 int digID = mDigits->size();
225 mDigits->emplace_back(chip.getChipIndex(), preDig.row, preDig.col, preDig.charge);
226 LOG(
debug) <<
"Adding digit ID: " << digID <<
" with chipID: " << chip.getChipIndex() <<
", row: " << preDig.row <<
", col: " << preDig.col <<
", charge: " << preDig.charge;
227 mMCLabels->
addElement(digID, preDig.labelRef.label);
228 auto& nextRef = preDig.labelRef;
229 while (nextRef.next >= 0) {
230 nextRef = extra[nextRef.next];
235 buffer.erase(itBeg, iter);
245 mROFRecords->push_back(rcROF);
249 mExtraBuff.emplace_back(mExtraBuff.front().release());
250 mExtraBuff.pop_front();
255void Digitizer::processHit(
const o2::itsmft::Hit& hit, uint32_t& maxFr,
int evID,
int srcID)
258 int subDetID = mGeometry->getSubDetID(chipID);
260 int layer = mGeometry->getLayer(chipID);
261 int disk = mGeometry->getDisk(chipID);
264 LOG(
debug) <<
"Skipping disk " << disk;
268 LOG(
debug) <<
"Processing hit for chip " << chipID;
269 auto& chip = mChips[chipID];
270 if (chip.isDisabled()) {
271 LOG(
debug) <<
"Skipping disabled chip " << chipID;
274 float timeInROF = hit.
GetTime() * sec2ns;
275 LOG(
debug) <<
"timeInROF: " << timeInROF;
276 if (timeInROF > 20e3) {
277 const int maxWarn = 10;
278 static int warnNo = 0;
279 if (warnNo < maxWarn) {
280 LOG(warning) <<
"Ignoring hit with time_in_event = " << timeInROF <<
" ns"
281 << ((++warnNo < maxWarn) ?
"" :
" (suppressing further warnings)");
286 timeInROF += mCollisionTimeWrtROF;
290 LOG(
debug) <<
"Ignoring hit with timeInROF = " << timeInROF;
303 int nFrames = roFrameRelMax + 1 - roFrameRel;
304 uint32_t roFrameMax = mNewROFrame + roFrameRelMax;
305 if (roFrameMax > maxFr) {
324 LOG(
debug) <<
"Called curved to flat: " << xyzLocS.x() <<
" -> " << xyFlatS.x() <<
", " << xyzLocS.y() <<
" -> " << xyFlatS.y();
326 xyzLocS.SetXYZ(xyFlatS.x(), xyFlatS.y(), xyzLocS.Z());
327 xyzLocE.SetXYZ(xyFlatE.x(), xyFlatE.y(), xyzLocE.Z());
337 LOG(
debug) <<
"local original: startPos = " << xyzLocS <<
", endPos = " << xyzLocE << std::endl;
338 xyzLocS.SetY(xyzLocS.Y());
339 xyzLocE.SetY(xyzLocE.Y());
341 LOG(
debug) <<
"rescaled Y: startPos = " << xyzLocS <<
", endPos = " << xyzLocE << std::endl;
352 LOG(
debug) <<
"Step into the sensitive volume: " <<
step <<
". Number of steps: " << nSteps;
353 int rowS = -1, colS = -1, rowE = -1, colE = -1, nSkip = 0;
358 if (++nSkip >= nSteps) {
359 LOG(
debug) <<
"Did not enter to sensitive matrix, " << nSkip <<
" >= " << nSteps;
367 if (++nSkip >= nSteps) {
368 LOG(
debug) <<
"Did not enter to sensitive matrix, " << nSkip <<
" >= " << nSteps;
379 std::swap(rowS, rowE);
382 std::swap(colS, colE);
400 int rowSpan = rowE - rowS + 1, colSpan = colE - colS + 1;
402 float respMatrix[rowSpan][colSpan];
403 std::fill(&respMatrix[0][0], &respMatrix[0][0] + rowSpan * colSpan, 0.f);
406 nElectrons *= nStepsInv;
411 int rowPrev = -1, colPrev = -1,
row,
col;
412 float cRowPix = 0.f, cColPix = 0.f;
421 xyzLocS.SetY(xyzLocS.Y() + ((subDetID == 0) ? mSimRespVDShift : mSimRespMLOTShift));
424 for (
int iStep = nSteps; iStep--;) {
427 if (
row != rowPrev ||
col != colPrev) {
434 bool flipCol =
false, flipRow =
false;
436 float rowMax{}, colMax{};
439 rowMax = 0.5f * Segmentation::PitchRowVD * mSimRespVDScaleX;
440 colMax = 0.5f * Segmentation::PitchColVD * mSimRespVDScaleZ;
441 rspmat = resp->
getResponse(mSimRespVDScaleX * (xyzLocS.X() - cRowPix), mSimRespVDScaleZ * (xyzLocS.Z() - cColPix), xyzLocS.Y(), flipRow, flipCol, rowMax, colMax);
443 rowMax = 0.5f * Segmentation::PitchRowMLOT * mSimRespMLOTScaleX;
444 colMax = 0.5f * Segmentation::PitchColMLOT * mSimRespMLOTScaleZ;
445 rspmat = resp->
getResponse(mSimRespMLOTScaleX * (xyzLocS.X() - cRowPix), mSimRespMLOTScaleZ * (xyzLocS.Z() - cColPix), xyzLocS.Y(), flipRow, flipCol, rowMax, colMax);
450 if (rspmat ==
nullptr) {
451 LOG(
debug) <<
"Error in rspmat for step " << iStep <<
" / " << nSteps;
454 LOG(
debug) <<
"rspmat valid! for step " << iStep <<
" / " << nSteps <<
", (row,col) = (" <<
row <<
"," <<
col <<
")";
459 if (rowDest < 0 || rowDest >= rowSpan) {
464 if (colDest < 0 || colDest >= colSpan) {
467 respMatrix[rowDest][colDest] += rspmat->getValue(irow, icol, mSimRespOrientation ? !flipRow : flipRow, !flipCol);
474 auto roFrameAbs = mNewROFrame + roFrameRel;
475 LOG(
debug) <<
"Spanning through rows and columns; rowspan = " << rowSpan <<
" colspan = " << colSpan <<
" = " << colE <<
" - " << colS <<
" +1 " << std::endl;
476 for (
int irow = rowSpan; irow--;) {
477 uint16_t rowIS = irow + rowS;
478 for (
int icol = colSpan; icol--;) {
479 float nEleResp = respMatrix[irow][icol];
480 if (nEleResp <= 1.e-36) {
483 LOG(
debug) <<
"nEleResp: value " << nEleResp <<
" for pixel " << irow <<
" " << icol << std::endl;
484 int nEle = gRandom->Poisson(nElectrons * nEleResp);
485 LOG(
debug) <<
"Charge detected in the pixel: " << nEle <<
" for pixel " << irow <<
" " << icol << std::endl;
488 LOG(
debug) <<
"Ignoring pixel with nEle = " << nEle <<
" < min charge to account "
493 uint16_t colIS = icol + colS;
494 if (mNoiseMap && mNoiseMap->
isNoisy(chipID, rowIS, colIS)) {
497 if (mDeadChanMap && mDeadChanMap->
isNoisy(chipID, rowIS, colIS)) {
501 registerDigits(chip, roFrameAbs, timeInROF, nFrames, rowIS, colIS, nEle, lbl);
513 LOG(
debug) <<
"Registering digits for chip " << chip.
getChipIndex() <<
" at ROFrame " << roFrame
514 <<
" row " <<
row <<
" col " <<
col <<
" nEle " << nEle <<
" label " << lbl;
516 for (
int i = 0;
i < nROF;
i++) {
517 uint32_t roFr = roFrame +
i;
525 if (roFr > mEventROFrameMax) {
526 mEventROFrameMax = roFr;
528 if (roFr < mEventROFrameMin) {
529 mEventROFrameMin = roFr;
535 LOG(
debug) <<
"Added digit " <<
key <<
" " << roFr <<
" " <<
row <<
" " <<
col <<
" " << nEleROF;
541 ExtraDig* extra = getExtraDigBuffer(roFr);
545 if ((*extra)[nxt].
label == lbl) {
549 nxt = (*extra)[nxt].next;
556 extra->emplace_back(lbl);
Definition of the SegmentationChipclass.
Definition of the TRK digitizer.
math_utils::Point3D< T > GetPos() const
short GetDetectorID() const
static const DPLDigitizerParam< N > & Instance()
const char * getName() const
const o2::detectors::DetID & getDetID() const
const Mat3D & getMatrixL2G(int sensID) const
static int constexpr NPix
float getMaxDuration() const
float getCollectedCharge(float totalNEle, float tMin, float tMax) const
bool getResponse(float vRow, float vCol, float cDepth, AlpideRespSimMat &dest) const
void addDigit(ULong64_t key, UInt_t roframe, UShort_t row, UShort_t col, int charge, o2::MCCompLabel lbl)
o2::itsmft::PreDigit * findDigit(ULong64_t key)
static ULong64_t getOrderingKey(UInt_t roframe, UShort_t row, UShort_t col)
Get global ordering key made of readout frame, column and row.
UShort_t getChipIndex() const
int getMinChargeToAccount() const
float getROFrameLengthInv() const
virtual void print() const
const SignalShape & getSignalShape() const
float getEnergyToNElectrons() const
int getROFrameLengthInBC() const
float getStrobeDelay() const
bool isContinuous() const
int getChargeThreshold() const
float getNSimStepsInv() const
float getROFrameLength() const
float getStrobeLength() const
void setEventTime(const o2::InteractionTimeRecord &irt)
auto getChipResponse(int chipID)
bool isContinuous() const
void process(const std::vector< Hit > *hits, int evID, int srcID)
Steer conversion of hits to digits.
void fillOutputContainer(uint32_t maxFrame=0xffffffff)
Int_t getNumberOfChips() const
math_utils::Point3D< Float_t > GetPosStart() const
bool isFullChipMasked(int chip) const
bool isNoisy(int chip, int row, int col) const
const BCData & getBCData() const
void setFirstEntry(int idx)
int getFirstEntry() const
void setROFrame(ROFtype rof)
static bool localToDetector(float x, float z, int &iRow, int &iCol)
static bool detectorToLocal(L row, L col, T &xRow, T &zCol)
static constexpr float PitchColVD
static constexpr float PitchColMLOT
static constexpr float PitchRowMLOT
static constexpr float PitchRowVD
GLuint GLsizei const GLchar * label
GLenum GLuint GLint GLint layer
constexpr double LHCBunchSpacingNS
constexpr std::array< int, nLayers > nRows
int64_t differenceInBC(const InteractionRecord &other) const
void setFromLong(int64_t l)
double timeInBCNS
time in NANOSECONDS relative to orbit/bc
int next
eventual next contribution to the same pixel
o2::MCCompLabel label
hit label
PreDigitLabelRef labelRef
label and reference to the next one
IR getFirstSampledTFIR() const
get TF and HB (abs) for this IR
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"