@@ -62,22 +62,29 @@ void Digitizer::init()
6262
6363 // / setting scale factors to adapt to the APTS response function (adjusting pitch and Y shift)
6464 // TODO: adjust Y shift when the geometry is improved
65+ LOG (debug) << " Depth max: " << mChipSimRespVD ->getDepthMax ();
66+ LOG (debug) << " Depth min: " << mChipSimRespVD ->getDepthMin ();
67+
68+ float thicknessVD = 0.0095 ; // cm --- hardcoded based on geometry currently present
69+ float thicknessMLOT = 0.1 ; // cm --- hardcoded based on geometry currently present
70+
6571 mSimRespVDScaleX = o2::trk::constants::apts::pitchX / o2::trk::SegmentationChip::PitchRowVD;
6672 mSimRespVDScaleZ = o2::trk::constants::apts::pitchZ / o2::trk::SegmentationChip::PitchColVD;
67- mSimRespVDScaleDepth = o2::trk::constants::apts::thickness / (0.1 ); // / introducing this scaling factor because the silicon thickness for the moment is 1 mm -> rescale to 45 um which is the depth of the APTS response
68- mSimRespVDShift = mChipSimRespVD ->getDepthMax () - 0.1 * mSimRespVDScaleDepth / 2 .f ; // the shift should be done considering the rescaling done to adapt to the wrong silicon thickness. TODO: remove the scaling factor for the depth when the silicon thickness match the simulated response
73+ mSimRespVDScaleDepth = o2::trk::constants::apts::thickness / (thicknessVD); // / introducing this scaling factor because the silicon thickness for the moment is 1 mm -> rescale to 45 um which is the depth of the APTS response
74+ // mSimRespVDShift = mChipSimRespVD->getDepthMax() - thicknessVD * mSimRespVDScaleDepth / 2.f; // the shift should be done considering the rescaling done to adapt to the wrong silicon thickness. TODO: remove the scaling factor for the depth when the silicon thickness match the simulated response
75+ mSimRespVDShift = mChipSimRespVD ->getDepthMax (); // the curved, rescaled, sensors have a width from 0 to -45. Must add 10 um (= max depth) to match the APTS response.
6976 mSimRespMLOTScaleX = o2::trk::constants::apts::pitchX / o2::trk::SegmentationChip::PitchRowMLOT;
7077 mSimRespMLOTScaleZ = o2::trk::constants::apts::pitchZ / o2::trk::SegmentationChip::PitchColMLOT;
71- mSimRespMLOTScaleDepth = o2::trk::constants::apts::thickness / (0.1 ); // / introducing this scaling factor because the silicon thickness for the moment is 1 mm -> rescale to 45 um which is the depth of the APTS response
72- mSimRespMLOTShift = mChipSimRespMLOT ->getDepthMax () - 0.1 * mSimRespMLOTScaleDepth / 2 .f ; // the shift should be done considering the rescaling done to adapt to the wrong silicon thickness. TODO: remove the scaling factor for the depth when the silicon thickness match the simulated response
73- mSimRespOrientation = true ;
78+ mSimRespMLOTScaleDepth = o2::trk::constants::apts::thickness / (thicknessMLOT ); // / introducing this scaling factor because the silicon thickness for the moment is 1 mm -> rescale to 45 um which is the depth of the APTS response
79+ mSimRespMLOTShift = mChipSimRespMLOT ->getDepthMax () - thicknessMLOT * mSimRespMLOTScaleDepth / 2 .f ; // the shift should be done considering the rescaling done to adapt to the wrong silicon thickness. TODO: remove the scaling factor for the depth when the silicon thickness match the simulated response
80+ mSimRespOrientation = false ;
7481
7582 // importing the parameters from DPLDigitizerParam.h
7683 auto & dOptTRK = DPLDigitizerParam<o2::detectors::DetID::TRK>::Instance ();
7784
7885 LOGP (info, " TRK Digitizer is initialised." );
7986 mParams .print ();
80- LOGP (info, " VD shift = {} = {} - {} ; ML/OT shift = {} = {} - {}" , mSimRespVDShift , mChipSimRespVD -> getDepthMax (), 0.1 * mSimRespVDScaleDepth / 2 . f , mSimRespMLOTShift , mChipSimRespMLOT ->getDepthMax (), 0.1 * mSimRespMLOTScaleDepth / 2 .f );
87+ LOGP (info, " VD shift = {} ; ML/OT shift = {} = {} - {}" , mSimRespVDShift , mSimRespMLOTShift , mChipSimRespMLOT ->getDepthMax (), thicknessMLOT * mSimRespMLOTScaleDepth / 2 .f );
8188 LOGP (info, " VD pixel scale on x = {} ; z = {}" , mSimRespVDScaleX , mSimRespVDScaleZ );
8289 LOGP (info, " ML/OT pixel scale on x = {} ; z = {}" , mSimRespMLOTScaleX , mSimRespMLOTScaleZ );
8390 LOGP (info, " Response orientation: {}" , mSimRespOrientation ? " flipped" : " normal" );
@@ -125,7 +132,6 @@ void Digitizer::process(const std::vector<Hit>* hits, int evID, int srcID)
125132 return (*hits)[lhs].GetDetectorID () < (*hits)[rhs].GetDetectorID ();
126133 });
127134 for (int i : hitIdx) {
128- // (*hits)[i].Print("");
129135 processHit ((*hits)[i], mROFrameMax , evID, srcID);
130136 }
131137
@@ -179,75 +185,78 @@ void Digitizer::setEventTime(const o2::InteractionTimeRecord& irt)
179185// _______________________________________________________________________
180186void Digitizer::fillOutputContainer (uint32_t frameLast)
181187{
182- // std::cout << "Entering fillOutputContainer " << std::endl;
183- // // // fill output with digits from min.cached up to requested frame, generating the noise beforehand
184- // if (frameLast > mROFrameMax) {
185- // frameLast = mROFrameMax;
186- // }
187- // // // make sure all buffers for extra digits are created up to the maxFrame
188- // getExtraDigBuffer(mROFrameMax);
189- // LOG(info) << "Filling " << mGeometry->getName() << " digits output for RO frames " << mROFrameMin << ":"
190- // << frameLast;
191-
192- // o2::itsmft::ROFRecord rcROF; /// using temporarly itsmft::ROFRecord
193-
194- // // we have to write chips in RO increasing order, therefore have to loop over the frames here
195- // for (; mROFrameMin <= frameLast; mROFrameMin++) {
196- // rcROF.setROFrame(mROFrameMin);
197- // rcROF.setFirstEntry(mDigits->size()); // start of current ROF in digits
198- // auto& extra = *(mExtraBuff.front().get());
199- // for (auto& chip : mChips) { /// TODO: do also for ML?
200- // if (chip.isDisabled()) {
201- // continue;
202- // }
203- // // chip.addNoiseVD(mROFrameMin, mROFrameMin, &mParams);
204- // auto& buffer = chip.getPreDigits();
205- // if (buffer.empty()) {
206- // continue;
207- // }
208- // auto itBeg = buffer.begin();
209- // auto iter = itBeg;
210- // ULong64_t maxKey = chip.getOrderingKey(mROFrameMin + 1, 0, 0) - 1; // fetch digits with key below that
211- // for (; iter != buffer.end(); ++iter) {
212- // if (iter->first > maxKey) {
213- // break; // is the digit ROFrame from the key > the max requested frame
214- // }
215- // auto& preDig = iter->second; // preDigit
216- // if (preDig.charge >= mParams.getChargeThreshold()) {
217- // int digID = mDigits->size();
218- // mDigits->emplace_back(chip.getChipIndex(), preDig.row, preDig.col, preDig.charge);
219- // mMCLabels->addElement(digID, preDig.labelRef.label);
220- // auto& nextRef = preDig.labelRef; // extra contributors are in extra array
221- // while (nextRef.next >= 0) {
222- // nextRef = extra[nextRef.next];
223- // mMCLabels->addElement(digID, nextRef.label);
224- // }
225- // }
226- // }
227- // buffer.erase(itBeg, iter);
228- // }
229- // // finalize ROF record
230- // rcROF.setNEntries(mDigits->size() - rcROF.getFirstEntry()); // number of digits
231- // if (isContinuous()) {
232- // rcROF.getBCData().setFromLong(mIRFirstSampledTF.toLong() + mROFrameMin * mParams.getROFrameLengthInBC());
233- // } else {
234- // rcROF.getBCData() = mEventTime; // RSTODO do we need to add trigger delay?
235- // }
236- // if (mROFRecords) {
237- // mROFRecords->push_back(rcROF);
238- // }
239- // extra.clear(); // clear container for extra digits of the mROFrameMin ROFrame
240- // // and move it as a new slot in the end
241- // mExtraBuff.emplace_back(mExtraBuff.front().release());
242- // mExtraBuff.pop_front();
243- // }
188+ std::cout << " Entering fillOutputContainer " << std::endl;
189+ std::cout << " Digit size before fill: " << mDigits ->size () << std::endl;
190+ // // fill output with digits from min.cached up to requested frame, generating the noise beforehand
191+ if (frameLast > mROFrameMax ) {
192+ frameLast = mROFrameMax ;
193+ }
194+ // // make sure all buffers for extra digits are created up to the maxFrame
195+ getExtraDigBuffer (mROFrameMax );
196+ LOG (info) << " Filling " << mGeometry ->getName () << " digits output for RO frames " << mROFrameMin << " :"
197+ << frameLast;
198+
199+ o2::itsmft::ROFRecord rcROF; // / using temporarly itsmft::ROFRecord
200+
201+ // we have to write chips in RO increasing order, therefore have to loop over the frames here
202+ for (; mROFrameMin <= frameLast; mROFrameMin ++) {
203+ std::cout << " Entering ROFrame " << mROFrameMin << std::endl;
204+ rcROF.setROFrame (mROFrameMin );
205+ rcROF.setFirstEntry (mDigits ->size ()); // start of current ROF in digits
206+
207+ auto & extra = *(mExtraBuff .front ().get ());
208+ for (auto & chip : mChips ) {
209+ if (chip.isDisabled ()) {
210+ continue ;
211+ }
212+ // chip.addNoise(mROFrameMin, mROFrameMin, &mParams); /// TODO: add noise
213+ auto & buffer = chip.getPreDigits ();
214+ if (buffer.empty ()) {
215+ continue ;
216+ }
217+ auto itBeg = buffer.begin ();
218+ auto iter = itBeg;
219+ ULong64_t maxKey = chip.getOrderingKey (mROFrameMin + 1 , 0 , 0 ) - 1 ; // fetch digits with key below that
220+ for (; iter != buffer.end (); ++iter) {
221+ if (iter->first > maxKey) {
222+ break ; // is the digit ROFrame from the key > the max requested frame
223+ }
224+ auto & preDig = iter->second ; // preDigit
225+ if (preDig.charge >= mParams .getChargeThreshold ()) {
226+ int digID = mDigits ->size ();
227+ mDigits ->emplace_back (chip.getChipIndex (), preDig.row , preDig.col , preDig.charge );
228+ LOG (debug) << " Adding digit ID: " << digID << " with chipID: " << chip.getChipIndex () << " , row: " << preDig.row << " , col: " << preDig.col << " , charge: " << preDig.charge ;
229+ mMCLabels ->addElement (digID, preDig.labelRef .label );
230+ auto & nextRef = preDig.labelRef ; // extra contributors are in extra array
231+ while (nextRef.next >= 0 ) {
232+ nextRef = extra[nextRef.next ];
233+ mMCLabels ->addElement (digID, nextRef.label );
234+ }
235+ }
236+ }
237+ buffer.erase (itBeg, iter);
238+ }
239+ // finalize ROF record
240+ rcROF.setNEntries (mDigits ->size () - rcROF.getFirstEntry ()); // number of digits
241+ if (isContinuous ()) {
242+ rcROF.getBCData ().setFromLong (mIRFirstSampledTF .toLong () + mROFrameMin * mParams .getROFrameLengthInBC ());
243+ } else {
244+ rcROF.getBCData () = mEventTime ; // RSTODO do we need to add trigger delay?
245+ }
246+ if (mROFRecords ->size ()) {
247+ mROFRecords ->push_back (rcROF);
248+ }
249+ extra.clear (); // clear container for extra digits of the mROFrameMin ROFrame
250+ // and move it as a new slot in the end
251+ mExtraBuff .emplace_back (mExtraBuff .front ().release ());
252+ mExtraBuff .pop_front ();
253+ }
244254}
245255
246256// _______________________________________________________________________
247257void Digitizer::processHit (const o2::itsmft::Hit& hit, uint32_t & maxFr, int evID, int srcID)
248258{
249259 int chipID = hit.GetDetectorID (); // // the chip ID at the moment is not referred to the chip but to a wider detector element (e.g. quarter of layer or disk in VD, stave in ML, half stave in OT)
250- LOG (debug) << " Processing hit for chip " << chipID << " in event " << evID << " from source " << srcID;
251260 int subDetID = mGeometry ->getSubDetID (chipID);
252261
253262 int layer = mGeometry ->getLayer (chipID);
@@ -257,6 +266,8 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID
257266 LOG (debug) << " Skipping disk " << disk;
258267 return ; // skipping hits on disks for the moment
259268 }
269+
270+ LOG (debug) << " Processing hit for chip " << chipID;
260271 auto & chip = mChips [chipID];
261272 if (chip.isDisabled ()) {
262273 LOG (debug) << " Skipping disabled chip " << chipID;
@@ -325,15 +336,15 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID
325336 // std::cout<<"Going back to glob coordinates: " << (matrix * exampleLoc) << std::endl;
326337
327338 // // adapting the depth (Y) of the chip to the APTS response maximum depth
339+ LOG (debug) << " local original: startPos = " << xyzLocS << " , endPos = " << xyzLocE << std::endl;
328340 if (subDetID == 0 ) {
329341 xyzLocS.SetY (xyzLocS.Y () * mSimRespVDScaleDepth );
330342 xyzLocE.SetY (xyzLocE.Y () * mSimRespVDScaleDepth );
331- LOG (debug) << " rescaled: startPos Y = " << xyzLocS.Y () << " , endPos Y = " << xyzLocE.Y () << std::endl;
332343 } else {
333344 xyzLocS.SetY (xyzLocS.Y () * mSimRespMLOTScaleDepth );
334345 xyzLocE.SetY (xyzLocE.Y () * mSimRespMLOTScaleDepth );
335- LOG (debug) << " rescaled: startPos Y = " << xyzLocS.Y () << " , endPos Y = " << xyzLocE.Y () << std::endl;
336346 }
347+ LOG (debug) << " rescaled Y: startPos = " << xyzLocS << " , endPos = " << xyzLocE << std::endl;
337348
338349 math_utils::Vector3D<float > step (xyzLocE);
339350 step -= xyzLocS;
@@ -344,7 +355,7 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID
344355 xyzLocS += stepH;
345356 xyzLocE -= stepH;
346357
347- LOG (debug) << " Step into the sensitive volume: " << step;
358+ LOG (debug) << " Step into the sensitive volume: " << step << " . Number of steps: " << nSteps ;
348359 int rowS = -1 , colS = -1 , rowE = -1 , colE = -1 , nSkip = 0 ;
349360
350361 // / here it is the control whether the hit is in the sensitive matrix based on the segmentation
@@ -366,8 +377,6 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID
366377 xyzLocE -= step;
367378 }
368379
369- LOG (info) << " rowS: " << rowS << " colS: " << colS << " rowE: " << rowE << " colE: " << colE;
370-
371380 int nCols = getNCols (subDetID, layer);
372381 int nRows = getNRows (subDetID, layer);
373382
@@ -428,7 +437,7 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID
428437 rowPrev = row;
429438 colPrev = col;
430439 }
431- bool flipCol = true , flipRow = true ;
440+ bool flipCol = false , flipRow = false ;
432441 // note that response needs coordinates along column row (locX) (locZ) then depth (locY)
433442 float rowMax{}, colMax{};
434443 const AlpideRespSimMat* rspmat{nullptr };
@@ -442,13 +451,24 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID
442451 rspmat = resp->getResponse (mSimRespMLOTScaleX * (xyzLocS.X () - cRowPix), mSimRespMLOTScaleZ * (xyzLocS.Z () - cColPix), xyzLocS.Y (), flipRow, flipCol, rowMax, colMax);
443452 }
444453
454+ float tempPitchX = 0 , tempPitchZ = 0 ;
455+ if (subDetID == 0 ) {
456+ tempPitchX = Segmentation::PitchRowVD;
457+ tempPitchZ = Segmentation::PitchColVD;
458+ } else {
459+ tempPitchX = Segmentation::PitchRowMLOT;
460+ tempPitchZ = Segmentation::PitchColMLOT;
461+ }
462+ LOG (debug) << " X and Z inside pixel at start = " << (xyzLocS.X () - cRowPix) << " , " << (xyzLocS.Z () - cColPix) << " , rescaled: " << mSimRespMLOTScaleX * (xyzLocS.X () - cRowPix) << " , " << mSimRespMLOTScaleZ * (xyzLocS.Z () - cColPix);
463+ LOG (debug) << " Hit inside pitch? X: " << ((xyzLocS.X () - cRowPix) < tempPitchX) << " Z: " << ((xyzLocS.Z () - cColPix) < tempPitchZ);
464+
445465 xyzLocS += step;
446466
447467 if (rspmat == nullptr ) {
448- LOG (error ) << " Error in rspmat for step " << iStep << " / " << nSteps;
468+ LOG (debug ) << " Error in rspmat for step " << iStep << " / " << nSteps;
449469 continue ;
450470 }
451- LOG (debug) << " rspmat valid! for step " << iStep << " / " << nSteps;
471+ LOG (debug) << " rspmat valid! for step " << iStep << " / " << nSteps << " , (row,col) = ( " << row << " , " << col << " ) " ;
452472 // rspmat->print(); // print the response matrix for debugging
453473
454474 for (int irow = AlpideRespSimMat::NPix; irow--;) {
@@ -469,9 +489,10 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID
469489 // fire the pixels assuming Poisson(n_response_electrons)
470490 o2::MCCompLabel lbl (hit.GetTrackID (), evID, srcID, false );
471491 auto roFrameAbs = mNewROFrame + roFrameRel;
472- for (int irow = rowSpan; irow--;) {
473- uint16_t rowIS = irow + rowS;
474- for (int icol = colSpan; icol--;) {
492+ LOG (debug) << " Spanning through rows and columns; rowspan = " << rowSpan << " colspan = " << colSpan << " = " << colE << " - " << colS << " +1 " << std::endl;
493+ for (int irow = rowSpan; irow--;) { // irow ranging from 4 to 0
494+ uint16_t rowIS = irow + rowS; // row distant irow from the row of the hit start
495+ for (int icol = colSpan; icol--;) { // icol ranging from 4 to 0
475496 float nEleResp = respMatrix[irow][icol]; // value of the probability of the response in this pixel
476497 if (nEleResp <= 1 .e -36 ) {
477498 continue ;
@@ -485,7 +506,8 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID
485506 << mParams .getMinChargeToAccount () << " for pixel " << irow << " " << icol;
486507 continue ;
487508 }
488- uint16_t colIS = icol + colS;
509+
510+ uint16_t colIS = icol + colS; // col distant icol from the col of the hit start
489511 if (mNoiseMap && mNoiseMap ->isNoisy (chipID, rowIS, colIS)) {
490512 continue ;
491513 }
0 commit comments