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(
debug) <<
" Depth max: " << mChipSimRespVD->getDepthMax();
66 LOG(
debug) <<
" Depth min: " << mChipSimRespVD->getDepthMin();
68 float thicknessVD = 0.0095;
69 float thicknessMLOT = 0.1;
75 mSimRespVDShift = mChipSimRespVD->getDepthMax();
79 mSimRespMLOTShift = mChipSimRespMLOT->getDepthMax() - thicknessMLOT * mSimRespMLOTScaleDepth / 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 * mSimRespMLOTScaleDepth / 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;
118 std::cout <<
"Printing segmentation info: " << std::endl;
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;
339 xyzLocS.SetY(xyzLocS.Y() * mSimRespVDScaleDepth);
340 xyzLocE.SetY(xyzLocE.Y() * mSimRespVDScaleDepth);
342 xyzLocS.SetY(xyzLocS.Y() * mSimRespMLOTScaleDepth);
343 xyzLocE.SetY(xyzLocE.Y() * mSimRespMLOTScaleDepth);
345 LOG(
debug) <<
"rescaled Y: startPos = " << xyzLocS <<
", endPos = " << xyzLocE << std::endl;
356 LOG(
debug) <<
"Step into the sensitive volume: " <<
step <<
". Number of steps: " << nSteps;
357 int rowS = -1, colS = -1, rowE = -1, colE = -1, nSkip = 0;
362 if (++nSkip >= nSteps) {
363 LOG(
debug) <<
"Did not enter to sensitive matrix, " << nSkip <<
" >= " << nSteps;
371 if (++nSkip >= nSteps) {
372 LOG(
debug) <<
"Did not enter to sensitive matrix, " << nSkip <<
" >= " << nSteps;
383 std::swap(rowS, rowE);
386 std::swap(colS, colE);
404 int rowSpan = rowE - rowS + 1, colSpan = colE - colS + 1;
406 float respMatrix[rowSpan][colSpan];
407 std::fill(&respMatrix[0][0], &respMatrix[0][0] + rowSpan * colSpan, 0.f);
410 nElectrons *= nStepsInv;
415 int rowPrev = -1, colPrev = -1,
row,
col;
416 float cRowPix = 0.f, cColPix = 0.f;
425 xyzLocS.SetY(xyzLocS.Y() + ((subDetID == 0) ? mSimRespVDShift : mSimRespMLOTShift));
428 for (
int iStep = nSteps; iStep--;) {
431 if (
row != rowPrev ||
col != colPrev) {
438 bool flipCol =
false, flipRow =
false;
440 float rowMax{}, colMax{};
443 rowMax = 0.5f * Segmentation::PitchRowVD * mSimRespVDScaleX;
444 colMax = 0.5f * Segmentation::PitchColVD * mSimRespVDScaleZ;
445 rspmat = resp->
getResponse(mSimRespVDScaleX * (xyzLocS.X() - cRowPix), mSimRespVDScaleZ * (xyzLocS.Z() - cColPix), xyzLocS.Y(), flipRow, flipCol, rowMax, colMax);
447 rowMax = 0.5f * Segmentation::PitchRowMLOT * mSimRespMLOTScaleX;
448 colMax = 0.5f * Segmentation::PitchColMLOT * mSimRespMLOTScaleZ;
449 rspmat = resp->
getResponse(mSimRespMLOTScaleX * (xyzLocS.X() - cRowPix), mSimRespMLOTScaleZ * (xyzLocS.Z() - cColPix), xyzLocS.Y(), flipRow, flipCol, rowMax, colMax);
452 float tempPitchX = 0, tempPitchZ = 0;
454 tempPitchX = Segmentation::PitchRowVD;
455 tempPitchZ = Segmentation::PitchColVD;
457 tempPitchX = Segmentation::PitchRowMLOT;
458 tempPitchZ = Segmentation::PitchColMLOT;
460 LOG(
debug) <<
"X and Z inside pixel at start = " << (xyzLocS.X() - cRowPix) <<
" , " << (xyzLocS.Z() - cColPix) <<
", rescaled: " << mSimRespMLOTScaleX * (xyzLocS.X() - cRowPix) <<
" , " << mSimRespMLOTScaleZ * (xyzLocS.Z() - cColPix);
461 LOG(
debug) <<
"Hit inside pitch? X: " << ((xyzLocS.X() - cRowPix) < tempPitchX) <<
" Z: " << ((xyzLocS.Z() - cColPix) < tempPitchZ);
465 if (rspmat ==
nullptr) {
466 LOG(
debug) <<
"Error in rspmat for step " << iStep <<
" / " << nSteps;
469 LOG(
debug) <<
"rspmat valid! for step " << iStep <<
" / " << nSteps <<
", (row,col) = (" <<
row <<
"," <<
col <<
")";
474 if (rowDest < 0 || rowDest >= rowSpan) {
479 if (colDest < 0 || colDest >= colSpan) {
482 respMatrix[rowDest][colDest] += rspmat->getValue(irow, icol, mSimRespOrientation ? !flipRow : flipRow, flipCol);
489 auto roFrameAbs = mNewROFrame + roFrameRel;
490 LOG(
debug) <<
"Spanning through rows and columns; rowspan = " << rowSpan <<
" colspan = " << colSpan <<
" = " << colE <<
" - " << colS <<
" +1 " << std::endl;
491 for (
int irow = rowSpan; irow--;) {
492 uint16_t rowIS = irow + rowS;
493 for (
int icol = colSpan; icol--;) {
494 float nEleResp = respMatrix[irow][icol];
495 if (nEleResp <= 1.e-36) {
498 LOG(
debug) <<
"nEleResp: value " << nEleResp <<
" for pixel " << irow <<
" " << icol << std::endl;
499 int nEle = gRandom->Poisson(nElectrons * nEleResp);
500 LOG(
debug) <<
"Charge detected in the pixel: " << nEle <<
" for pixel " << irow <<
" " << icol << std::endl;
503 LOG(
debug) <<
"Ignoring pixel with nEle = " << nEle <<
" < min charge to account "
508 uint16_t colIS = icol + colS;
509 if (mNoiseMap && mNoiseMap->
isNoisy(chipID, rowIS, colIS)) {
512 if (mDeadChanMap && mDeadChanMap->
isNoisy(chipID, rowIS, colIS)) {
516 registerDigits(chip, roFrameAbs, timeInROF, nFrames, rowIS, colIS, nEle, lbl);
528 LOG(
debug) <<
"Registering digits for chip " << chip.
getChipIndex() <<
" at ROFrame " << roFrame
529 <<
" row " <<
row <<
" col " <<
col <<
" nEle " << nEle <<
" label " << lbl;
531 for (
int i = 0;
i < nROF;
i++) {
532 uint32_t roFr = roFrame +
i;
540 if (roFr > mEventROFrameMax) {
541 mEventROFrameMax = roFr;
543 if (roFr < mEventROFrameMin) {
544 mEventROFrameMin = roFr;
550 LOG(
debug) <<
"Added digit " <<
key <<
" " << roFr <<
" " <<
row <<
" " <<
col <<
" " << nEleROF;
556 ExtraDig* extra = getExtraDigBuffer(roFr);
560 if ((*extra)[nxt].
label == lbl) {
564 nxt = (*extra)[nxt].next;
571 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 const void Print() noexcept
Print segmentation info.
static constexpr float PitchRowVD
GLuint GLsizei const GLchar * label
GLenum GLuint GLint GLint layer
constexpr double LHCBunchSpacingNS
constexpr std::array< int, nLayers > nRows
constexpr double thickness
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"