diff --git a/Detectors/Upgrades/ALICE3/TRK/CMakeLists.txt b/Detectors/Upgrades/ALICE3/TRK/CMakeLists.txt index 83838a01d13f1..645e3149e4ab7 100644 --- a/Detectors/Upgrades/ALICE3/TRK/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/TRK/CMakeLists.txt @@ -10,4 +10,5 @@ # or submit itself to any jurisdiction. add_subdirectory(base) -add_subdirectory(simulation) \ No newline at end of file +add_subdirectory(simulation) +add_subdirectory(workflow) \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/CMakeLists.txt b/Detectors/Upgrades/ALICE3/TRK/simulation/CMakeLists.txt index 7706c0e10d778..9a775fa87ce57 100644 --- a/Detectors/Upgrades/ALICE3/TRK/simulation/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/CMakeLists.txt @@ -12,12 +12,16 @@ o2_add_library(TRKSimulation SOURCES src/TRKLayer.cxx src/Detector.cxx + src/Digitizer.cxx src/TRKServices.cxx + src/DPLDigitizerParam.cxx PUBLIC_LINK_LIBRARIES O2::TRKBase O2::FT3Simulation - O2::ITSMFTSimulation) + O2::ITSMFTSimulation + O2::SimulationDataFormat) o2_target_root_dictionary(TRKSimulation - HEADERS include/TRKSimulation/Detector.h + HEADERS include/TRKSimulation/Digitizer.h + include/TRKSimulation/Detector.h include/TRKSimulation/TRKLayer.h include/TRKSimulation/TRKServices.h) \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DPLDigitizerParam.h b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DPLDigitizerParam.h new file mode 100644 index 0000000000000..59b3551ecbd32 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/DPLDigitizerParam.h @@ -0,0 +1,69 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef ALICEO2_TRKDPLDIGITIZERPARAM_H_ +#define ALICEO2_TRKDPLDIGITIZERPARAM_H_ + +#include "DetectorsCommonDataFormats/DetID.h" +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/ConfigurableParamHelper.h" +#include + +namespace o2 +{ +namespace trk +{ +template +struct DPLDigitizerParam : public o2::conf::ConfigurableParamHelper> { + static_assert(N == o2::detectors::DetID::TRK || N == o2::detectors::DetID::FT3, "only DetID::TRK or DetID::FT3 are allowed"); + + static constexpr std::string_view getParamName() + { + return N == o2::detectors::DetID::TRK ? ParamName[0] : ParamName[1]; + } + + bool continuous = true; ///< flag for continuous simulation + float noisePerPixel = DEFNoisePerPixel(); ///< ALPIDE Noise per channel + float strobeFlatTop = 7500.; ///< strobe shape flat top + float strobeMaxRiseTime = 1100.; ///< strobe max rise time + float strobeQRiseTime0 = 450.; ///< q @ which strobe rise time is 0 + + double timeOffset = 0.; ///< time offset (in seconds!) to calculate ROFrame from hit time + int chargeThreshold = 150; ///< charge threshold in Nelectrons + int minChargeToAccount = 15; ///< minimum charge contribution to account + int nSimSteps = 7; ///< number of steps in response simulation + float energyToNElectrons = 1. / 3.6e-9; // conversion of eloss to Nelectrons + + float Vbb = 0.0; ///< back bias absolute value for MFT (in Volt) + float IBVbb = 0.0; ///< back bias absolute value for ITS Inner Barrel (in Volt) + float OBVbb = 0.0; ///< back bias absolute value for ITS Outter Barrel (in Volt) + + std::string noiseFilePath{}; ///< optional noise masks file path. FIXME to be removed once switch to CCDBFetcher + + // boilerplate stuff + make principal key + O2ParamDef(DPLDigitizerParam, getParamName().data()); + + private: + static constexpr float DEFNoisePerPixel() + { + return N == o2::detectors::DetID::TRK ? 1e-8 : 1e-8; // ITS/MFT values here!! + } + + static constexpr std::string_view ParamName[2] = {"TRKDigitizerParam", "FT3DigitizerParam"}; +}; + +template +DPLDigitizerParam DPLDigitizerParam::sInstance; + +} // namespace trk +} // namespace o2 + +#endif diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/Digitizer.h b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/Digitizer.h new file mode 100644 index 0000000000000..6863c5392cae3 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/Digitizer.h @@ -0,0 +1,128 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file Digitizer.h +/// \brief Definition of the TRK digitizer +#ifndef ALICEO2_TRK_DIGITIZER_H +#define ALICEO2_TRK_DIGITIZER_H + +#include +#include +#include + +#include "Rtypes.h" // for Digitizer::Class +#include "TObject.h" // for TObject + +#include "ITSMFTSimulation/ChipDigitsContainer.h" +// #include "ITSMFTSimulation/AlpideSimResponse.h" +#include "ITSMFTSimulation/DigiParams.h" +#include "ITSMFTSimulation/Hit.h" +#include "TRKBase/GeometryTGeo.h" +// #include "ITS3Base/SegmentationSuperAlpide.h" +#include "DataFormatsITSMFT/Digit.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "CommonDataFormat/InteractionRecord.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#endif + +namespace o2::trk +{ + +class Digitizer : public TObject +{ + using ExtraDig = std::vector; ///< container for extra contributions to PreDigits + + public: + void setDigits(std::vector* dig) { mDigits = dig; } + void setMCLabels(o2::dataformats::MCTruthContainer* mclb) { mMCLabels = mclb; } + void setROFRecords(std::vector* rec) { mROFRecords = rec; } + + o2::itsmft::DigiParams& getParams() { return (o2::itsmft::DigiParams&)mParams; } + const o2::itsmft::DigiParams& getParams() const { return mParams; } + + void init(); + + /// Steer conversion of hits to digits + void process(const std::vector* hits, int evID, int srcID); + void setEventTime(const o2::InteractionTimeRecord& irt); + double getEndTimeOfROFMax() const + { + ///< return the time corresponding to end of the last reserved ROFrame : mROFrameMax + return mParams.getROFrameLength() * (mROFrameMax + 1) + mParams.getTimeOffset(); + } + + void setContinuous(bool v) { mParams.setContinuous(v); } + bool isContinuous() const { return mParams.isContinuous(); } + void fillOutputContainer(uint32_t maxFrame = 0xffffffff); + + void setDigiParams(const o2::itsmft::DigiParams& par) { mParams = par; } + const o2::itsmft::DigiParams& getDigitParams() const { return mParams; } + + // provide the common itsmft::GeometryTGeo to access matrices and segmentation + void setGeometry(const o2::trk::GeometryTGeo* gm) { mGeometry = gm; } + + uint32_t getEventROFrameMin() const { return mEventROFrameMin; } + uint32_t getEventROFrameMax() const { return mEventROFrameMax; } + void resetEventROFrames() + { + mEventROFrameMin = 0xffffffff; + mEventROFrameMax = 0; + } + + void setDeadChannelsMap(const o2::itsmft::NoiseMap* mp) { mDeadChanMap = mp; } + + private: + void processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID, int srcID); + void registerDigits(o2::itsmft::ChipDigitsContainer& chip, uint32_t roFrame, float tInROF, int nROF, + uint16_t row, uint16_t col, int nEle, o2::MCCompLabel& lbl); + + ExtraDig* getExtraDigBuffer(uint32_t roFrame) + { + if (mROFrameMin > roFrame) { + return nullptr; // nothing to do + } + int ind = roFrame - mROFrameMin; + while (ind >= int(mExtraBuff.size())) { + mExtraBuff.emplace_back(std::make_unique()); + } + return mExtraBuff[ind].get(); + } + + static constexpr float sec2ns = 1e9; + + o2::itsmft::DigiParams mParams; ///< digitization parameters + o2::InteractionTimeRecord mEventTime; ///< global event time and interaction record + o2::InteractionRecord mIRFirstSampledTF; ///< IR of the 1st sampled IR, noise-only ROFs will be inserted till this IR only + double mCollisionTimeWrtROF{}; + uint32_t mROFrameMin = 0; ///< lowest RO frame of current digits + uint32_t mROFrameMax = 0; ///< highest RO frame of current digits + uint32_t mNewROFrame = 0; ///< ROFrame corresponding to provided time + + uint32_t mEventROFrameMin = 0xffffffff; ///< lowest RO frame for processed events (w/o automatic noise ROFs) + uint32_t mEventROFrameMax = 0; ///< highest RO frame forfor processed events (w/o automatic noise ROFs) + + o2::itsmft::AlpideSimResponse* mAlpSimResp = nullptr; // simulated response + + const o2::trk::GeometryTGeo* mGeometry = nullptr; ///< TRK geometry + + std::vector mChips; ///< Array of chips digits containers + std::deque> mExtraBuff; ///< burrer (per roFrame) for extra digits + + std::vector* mDigits = nullptr; //! output digits + std::vector* mROFRecords = nullptr; //! output ROF records + o2::dataformats::MCTruthContainer* mMCLabels = nullptr; //! output labels + + const o2::itsmft::NoiseMap* mDeadChanMap = nullptr; + + ClassDef(Digitizer, 1); +}; +} // namespace o2::trk \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/src/DPLDigitizerParam.cxx b/Detectors/Upgrades/ALICE3/TRK/simulation/src/DPLDigitizerParam.cxx new file mode 100644 index 0000000000000..a13f2e58bd3a4 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/src/DPLDigitizerParam.cxx @@ -0,0 +1,23 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "TRKSimulation/DPLDigitizerParam.h" + +namespace o2 +{ +namespace trk +{ +// this makes sure that the constructor of the parameters is statically called +// so that these params are part of the parameter database +static auto& sDigitizerParamITS = o2::trk::DPLDigitizerParam::Instance(); +static auto& sDigitizerParamMFT = o2::trk::DPLDigitizerParam::Instance(); +} // namespace trk +} // namespace o2 diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ALICE3/TRK/simulation/src/Digitizer.cxx new file mode 100644 index 0000000000000..21e6e629ec418 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/src/Digitizer.cxx @@ -0,0 +1,467 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file Digitizer.cxx + +#include "DataFormatsITSMFT/Digit.h" +// #include "ITSMFTBase/SegmentationAlpide.h" +#include "TRKSimulation/DPLDigitizerParam.h" +#include "TRKSimulation/Digitizer.h" +// #include "MathUtils/Cartesian.h" +// #include "SimulationDataFormat/MCTruthContainer.h" +// #include "DetectorsRaw/HBFUtils.h" + +// #include +// #include +// #include +// #include +#include // for LOG + +using o2::itsmft::Digit; +using o2::itsmft::Hit; +// using Segmentation = o2::itsmft::SegmentationAlpide; + +using namespace o2::trk; +// using namespace o2::base; + +//_______________________________________________________________________ +void Digitizer::init() +{ + // mNumberOfChips = mGeometry->getNumberOfChips(); + // mChips.resize(mNumberOfChips); + // for (int i = mNumberOfChips; i--;) { + // mChips[i].setChipIndex(i); + // if (mNoiseMap) { + // mChips[i].setNoiseMap(mNoiseMap); + // } + // if (mDeadChanMap) { + // mChips[i].disable(mDeadChanMap->isFullChipMasked(i)); + // mChips[i].setDeadChanMap(mDeadChanMap); + // } + // } + // initializing for both collection tables + /*for (int i = 0; i < 2; i++) { + mAlpSimResp[i].initData(i); + }*/ + + // importing the charge collection tables + // (initialized while building O2) + // auto file = TFile::Open(mResponseFile.data()); + // if (!file) { + // LOG(fatal) << "Cannot open response file " << mResponseFile; + // } + /*std::string response = "response"; + for (int i=0; i<2; i++) { + response.append(std::to_string(i)); + mAlpSimResp[i] = *(o2::itsmft::AlpideSimResponse*)file->Get(response.data()); + }*/ + // mAlpSimResp[0] = *(o2::itsmft::AlpideSimResponse*)file->Get("response0"); + // mAlpSimResp[1] = *(o2::itsmft::AlpideSimResponse*)file->Get("response1"); + + // importing the parameters from DPLDigitizerParam.h + auto& dOptTRK = DPLDigitizerParam::Instance(); + + LOGP(info, "TRK Digitizer is initalised."); +} + +// auto Digitizer::getChipResponse(int chipID) +// { +// if (mNumberOfChips < 10000) { // in MFT +// return mAlpSimRespMFT; +// } + +// if (chipID < 432) { // in ITS Inner Barrel +// return mAlpSimRespIB; +// } else { // in ITS Outter Barrel +// return mAlpSimRespOB; +// } +// } + +//_______________________________________________________________________ +void Digitizer::process(const std::vector* hits, int evID, int srcID) +{ + // digitize single event, the time must have been set beforehand + + // LOG(info) << "Digitizing " << mGeometry->getName() << " hits of entry " << evID << " from source " + // << srcID << " at time " << mEventTime << " ROFrame= " << mNewROFrame << ")" + // << " cont.mode: " << isContinuous() + // << " Min/Max ROFrames " << mROFrameMin << "/" << mROFrameMax; + + // // is there something to flush ? + // if (mNewROFrame > mROFrameMin) { + // fillOutputContainer(mNewROFrame - 1); // flush out all frame preceding the new one + // } + + // int nHits = hits->size(); + // std::vector hitIdx(nHits); + // std::iota(std::begin(hitIdx), std::end(hitIdx), 0); + // // sort hits to improve memory access + // std::sort(hitIdx.begin(), hitIdx.end(), + // [hits](auto lhs, auto rhs) { + // return (*hits)[lhs].GetDetectorID() < (*hits)[rhs].GetDetectorID(); + // }); + // for (int i : hitIdx) { + // processHit((*hits)[i], mROFrameMax, evID, srcID); + // } + // // in the triggered mode store digits after every MC event + // // TODO: in the real triggered mode this will not be needed, this is actually for the + // // single event processing only + // if (!mParams.isContinuous()) { + // fillOutputContainer(mROFrameMax); + // } +} + +//_______________________________________________________________________ +void Digitizer::setEventTime(const o2::InteractionTimeRecord& irt) +{ + // // assign event time in ns + // mEventTime = irt; + // if (!mParams.isContinuous()) { + // mROFrameMin = 0; // in triggered mode reset the frame counters + // mROFrameMax = 0; + // } + // // RO frame corresponding to provided time + // mCollisionTimeWrtROF = mEventTime.timeInBCNS; // in triggered mode the ROF starts at BC (is there a delay?) + // if (mParams.isContinuous()) { + // auto nbc = mEventTime.differenceInBC(mIRFirstSampledTF); + // if (mCollisionTimeWrtROF < 0 && nbc > 0) { + // nbc--; + // } + + // // we might get interactions to digitize from before + // // the first sampled IR + // if (nbc < 0) { + // mNewROFrame = 0; + // // this event is before the first RO + // mIsBeforeFirstRO = true; + // } else { + // mNewROFrame = nbc / mParams.getROFrameLengthInBC(); + // mIsBeforeFirstRO = false; + // } + // LOG(info) << " NewROFrame " << mNewROFrame << " nbc " << nbc; + + // // in continuous mode depends on starts of periodic readout frame + // mCollisionTimeWrtROF += (nbc % mParams.getROFrameLengthInBC()) * o2::constants::lhc::LHCBunchSpacingNS; + // } else { + // mNewROFrame = 0; + // } + + // if (mNewROFrame < mROFrameMin) { + // LOG(error) << "New ROFrame " << mNewROFrame << " (" << irt << ") precedes currently cashed " << mROFrameMin; + // throw std::runtime_error("deduced ROFrame precedes already processed one"); + // } + + // if (mParams.isContinuous() && mROFrameMax < mNewROFrame) { + // mROFrameMax = mNewROFrame - 1; // all frames up to this are finished + // } +} + +//_______________________________________________________________________ +void Digitizer::fillOutputContainer(uint32_t frameLast) +{ + // // fill output with digits from min.cached up to requested frame, generating the noise beforehand + // if (frameLast > mROFrameMax) { + // frameLast = mROFrameMax; + // } + // // make sure all buffers for extra digits are created up to the maxFrame + // getExtraDigBuffer(mROFrameMax); + + // LOG(info) << "Filling " << mGeometry->getName() << " digits output for RO frames " << mROFrameMin << ":" + // << frameLast; + + // o2::itsmft::ROFRecord rcROF; + + // // we have to write chips in RO increasing order, therefore have to loop over the frames here + // for (; mROFrameMin <= frameLast; mROFrameMin++) { + // rcROF.setROFrame(mROFrameMin); + // rcROF.setFirstEntry(mDigits->size()); // start of current ROF in digits + + // auto& extra = *(mExtraBuff.front().get()); + // for (auto& chip : mChips) { + // if (chip.isDisabled()) { + // continue; + // } + // chip.addNoise(mROFrameMin, mROFrameMin, &mParams); + // auto& buffer = chip.getPreDigits(); + // if (buffer.empty()) { + // continue; + // } + // auto itBeg = buffer.begin(); + // auto iter = itBeg; + // ULong64_t maxKey = chip.getOrderingKey(mROFrameMin + 1, 0, 0) - 1; // fetch digits with key below that + // for (; iter != buffer.end(); ++iter) { + // if (iter->first > maxKey) { + // break; // is the digit ROFrame from the key > the max requested frame + // } + // auto& preDig = iter->second; // preDigit + // if (preDig.charge >= mParams.getChargeThreshold()) { + // int digID = mDigits->size(); + // mDigits->emplace_back(chip.getChipIndex(), preDig.row, preDig.col, preDig.charge); + // mMCLabels->addElement(digID, preDig.labelRef.label); + // auto& nextRef = preDig.labelRef; // extra contributors are in extra array + // while (nextRef.next >= 0) { + // nextRef = extra[nextRef.next]; + // mMCLabels->addElement(digID, nextRef.label); + // } + // } + // } + // buffer.erase(itBeg, iter); + // } + // // finalize ROF record + // rcROF.setNEntries(mDigits->size() - rcROF.getFirstEntry()); // number of digits + // if (isContinuous()) { + // rcROF.getBCData().setFromLong(mIRFirstSampledTF.toLong() + mROFrameMin * mParams.getROFrameLengthInBC()); + // } else { + // rcROF.getBCData() = mEventTime; // RSTODO do we need to add trigger delay? + // } + // if (mROFRecords) { + // mROFRecords->push_back(rcROF); + // } + // extra.clear(); // clear container for extra digits of the mROFrameMin ROFrame + // // and move it as a new slot in the end + // mExtraBuff.emplace_back(mExtraBuff.front().release()); + // mExtraBuff.pop_front(); + // } +} + +//_______________________________________________________________________ +void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID, int srcID) +{ + // // convert single hit to digits + // int chipID = hit.GetDetectorID(); + // auto& chip = mChips[chipID]; + // if (chip.isDisabled()) { + // LOG(debug) << "skip disabled chip " << chipID; + // return; + // } + // float timeInROF = hit.GetTime() * sec2ns; + // if (timeInROF > 20e3) { + // const int maxWarn = 10; + // static int warnNo = 0; + // if (warnNo < maxWarn) { + // LOG(warning) << "Ignoring hit with time_in_event = " << timeInROF << " ns" + // << ((++warnNo < maxWarn) ? "" : " (suppressing further warnings)"); + // } + // return; + // } + // if (isContinuous()) { + // timeInROF += mCollisionTimeWrtROF; + // } + // if (mIsBeforeFirstRO && timeInROF < 0) { + // // disregard this hit because it comes from an event before readout starts and it does not effect this RO + // return; + // } + + // // calculate RO Frame for this hit + // if (timeInROF < 0) { + // timeInROF = 0.; + // } + // float tTot = mParams.getSignalShape().getMaxDuration(); + // // frame of the hit signal start wrt event ROFrame + // int roFrameRel = int(timeInROF * mParams.getROFrameLengthInv()); + // // frame of the hit signal end wrt event ROFrame: in the triggered mode we read just 1 frame + // uint32_t roFrameRelMax = mParams.isContinuous() ? (timeInROF + tTot) * mParams.getROFrameLengthInv() : roFrameRel; + // int nFrames = roFrameRelMax + 1 - roFrameRel; + // uint32_t roFrameMax = mNewROFrame + roFrameRelMax; + // if (roFrameMax > maxFr) { + // maxFr = roFrameMax; // if signal extends beyond current maxFrame, increase the latter + // } + + // // here we start stepping in the depth of the sensor to generate charge diffusion + // float nStepsInv = mParams.getNSimStepsInv(); + // int nSteps = mParams.getNSimSteps(); + // const auto& matrix = mGeometry->getMatrixL2G(hit.GetDetectorID()); + // math_utils::Vector3D xyzLocS(matrix ^ (hit.GetPosStart())); // start position in sensor frame + // math_utils::Vector3D xyzLocE(matrix ^ (hit.GetPos())); // end position in sensor frame + + // math_utils::Vector3D step(xyzLocE); + // step -= xyzLocS; + // step *= nStepsInv; // position increment at each step + // // the electrons will injected in the middle of each step + // math_utils::Vector3D stepH(step * 0.5); + // xyzLocS += stepH; + // xyzLocE -= stepH; + + // int rowS = -1, colS = -1, rowE = -1, colE = -1, nSkip = 0; + // // get entrance pixel row and col + // while (!Segmentation::localToDetector(xyzLocS.X(), xyzLocS.Z(), rowS, colS)) { // guard-ring ? + // if (++nSkip >= nSteps) { + // return; // did not enter to sensitive matrix + // } + // xyzLocS += step; + // } + // // get exit pixel row and col + // while (!Segmentation::localToDetector(xyzLocE.X(), xyzLocE.Z(), rowE, colE)) { // guard-ring ? + // if (++nSkip >= nSteps) { + // return; // did not enter to sensitive matrix + // } + // xyzLocE -= step; + // } + // // estimate the limiting min/max row and col where the non-0 response is possible + // if (rowS > rowE) { + // std::swap(rowS, rowE); + // } + // if (colS > colE) { + // std::swap(colS, colE); + // } + // rowS -= AlpideRespSimMat::NPix / 2; + // rowE += AlpideRespSimMat::NPix / 2; + // if (rowS < 0) { + // rowS = 0; + // } + // if (rowE >= Segmentation::NRows) { + // rowE = Segmentation::NRows - 1; + // } + // colS -= AlpideRespSimMat::NPix / 2; + // colE += AlpideRespSimMat::NPix / 2; + // if (colS < 0) { + // colS = 0; + // } + // if (colE >= Segmentation::NCols) { + // colE = Segmentation::NCols - 1; + // } + // int rowSpan = rowE - rowS + 1, colSpan = colE - colS + 1; // size of plaquet where some response is expected + + // float respMatrix[rowSpan][colSpan]; // response accumulated here + // std::fill(&respMatrix[0][0], &respMatrix[0][0] + rowSpan * colSpan, 0.f); + + // float nElectrons = hit.GetEnergyLoss() * mParams.getEnergyToNElectrons(); // total number of deposited electrons + // nElectrons *= nStepsInv; // N electrons injected per step + // if (nSkip) { + // nSteps -= nSkip; + // } + // // + // int rowPrev = -1, colPrev = -1, row, col; + // float cRowPix = 0.f, cColPix = 0.f; // local coordinated of the current pixel center + + // const o2::itsmft::AlpideSimResponse* resp = getChipResponse(chipID); + + // // take into account that the AlpideSimResponse depth defintion has different min/max boundaries + // // although the max should coincide with the surface of the epitaxial layer, which in the chip + // // local coordinates has Y = +SensorLayerThickness/2 + + // xyzLocS.SetY(xyzLocS.Y() + resp->getDepthMax() - Segmentation::SensorLayerThickness / 2.); + + // // collect charge in every pixel which might be affected by the hit + // for (int iStep = nSteps; iStep--;) { + // // Get the pixel ID + // Segmentation::localToDetector(xyzLocS.X(), xyzLocS.Z(), row, col); + // if (row != rowPrev || col != colPrev) { // update pixel and coordinates of its center + // if (!Segmentation::detectorToLocal(row, col, cRowPix, cColPix)) { + // continue; // should not happen + // } + // rowPrev = row; + // colPrev = col; + // } + // bool flipCol, flipRow; + // // note that response needs coordinates along column row (locX) (locZ) then depth (locY) + // auto rspmat = resp->getResponse(xyzLocS.X() - cRowPix, xyzLocS.Z() - cColPix, xyzLocS.Y(), flipRow, flipCol); + + // xyzLocS += step; + // if (!rspmat) { + // continue; + // } + + // for (int irow = AlpideRespSimMat::NPix; irow--;) { + // int rowDest = row + irow - AlpideRespSimMat::NPix / 2 - rowS; // destination row in the respMatrix + // if (rowDest < 0 || rowDest >= rowSpan) { + // continue; + // } + // for (int icol = AlpideRespSimMat::NPix; icol--;) { + // int colDest = col + icol - AlpideRespSimMat::NPix / 2 - colS; // destination column in the respMatrix + // if (colDest < 0 || colDest >= colSpan) { + // continue; + // } + // respMatrix[rowDest][colDest] += rspmat->getValue(irow, icol, flipRow, flipCol); + // } + // } + // } + + // // fire the pixels assuming Poisson(n_response_electrons) + // o2::MCCompLabel lbl(hit.GetTrackID(), evID, srcID, false); + // auto roFrameAbs = mNewROFrame + roFrameRel; + // for (int irow = rowSpan; irow--;) { + // uint16_t rowIS = irow + rowS; + // for (int icol = colSpan; icol--;) { + // float nEleResp = respMatrix[irow][icol]; + // if (!nEleResp) { + // continue; + // } + // int nEle = gRandom->Poisson(nElectrons * nEleResp); // total charge in given pixel + // // ignore charge which have no chance to fire the pixel + // if (nEle < mParams.getMinChargeToAccount()) { + // continue; + // } + // uint16_t colIS = icol + colS; + // if (mNoiseMap && mNoiseMap->isNoisy(chipID, rowIS, colIS)) { + // continue; + // } + // if (mDeadChanMap && mDeadChanMap->isNoisy(chipID, rowIS, colIS)) { + // continue; + // } + // // + // registerDigits(chip, roFrameAbs, timeInROF, nFrames, rowIS, colIS, nEle, lbl); + // } + // } +} + +//________________________________________________________________________________ +void Digitizer::registerDigits(o2::itsmft::ChipDigitsContainer& chip, uint32_t roFrame, float tInROF, int nROF, + uint16_t row, uint16_t col, int nEle, o2::MCCompLabel& lbl) +{ + // Register digits for given pixel, accounting for the possible signal contribution to + // multiple ROFrame. The signal starts at time tInROF wrt the start of provided roFrame + // In every ROFrame we check the collected signal during strobe + + // float tStrobe = mParams.getStrobeDelay() - tInROF; // strobe start wrt signal start + // for (int i = 0; i < nROF; i++) { + // uint32_t roFr = roFrame + i; + // int nEleROF = mParams.getSignalShape().getCollectedCharge(nEle, tStrobe, tStrobe + mParams.getStrobeLength()); + // tStrobe += mParams.getROFrameLength(); // for the next ROF + + // // discard too small contributions, they have no chance to produce a digit + // if (nEleROF < mParams.getMinChargeToAccount()) { + // continue; + // } + // if (roFr > mEventROFrameMax) { + // mEventROFrameMax = roFr; + // } + // if (roFr < mEventROFrameMin) { + // mEventROFrameMin = roFr; + // } + // auto key = chip.getOrderingKey(roFr, row, col); + // PreDigit* pd = chip.findDigit(key); + // if (!pd) { + // chip.addDigit(key, roFr, row, col, nEleROF, lbl); + // } else { // there is already a digit at this slot, account as PreDigitExtra contribution + // pd->charge += nEleROF; + // if (pd->labelRef.label == lbl) { // don't store the same label twice + // continue; + // } + // ExtraDig* extra = getExtraDigBuffer(roFr); + // int& nxt = pd->labelRef.next; + // bool skip = false; + // while (nxt >= 0) { + // if ((*extra)[nxt].label == lbl) { // don't store the same label twice + // skip = true; + // break; + // } + // nxt = (*extra)[nxt].next; + // } + // if (skip) { + // continue; + // } + // // new predigit will be added in the end of the chain + // nxt = extra->size(); + // extra->emplace_back(lbl); + // } + // } +} diff --git a/Detectors/Upgrades/ALICE3/TRK/simulation/src/TRKSimulationLinkDef.h b/Detectors/Upgrades/ALICE3/TRK/simulation/src/TRKSimulationLinkDef.h index 1b0181144b5d4..3f1da6257c3a0 100644 --- a/Detectors/Upgrades/ALICE3/TRK/simulation/src/TRKSimulationLinkDef.h +++ b/Detectors/Upgrades/ALICE3/TRK/simulation/src/TRKSimulationLinkDef.h @@ -19,5 +19,11 @@ #pragma link C++ class o2::trk::TRKServices + ; #pragma link C++ class o2::trk::Detector + ; #pragma link C++ class o2::base::DetImpl < o2::trk::Detector> + ; +#pragma link C++ class o2::trk::Digitizer + ; + +// #pragma link C++ class o2::itsmft::DPLDigitizerParam < o2::detectors::DetID::ITS> + ; +// #pragma link C++ class o2::itsmft::DPLDigitizerParam < o2::detectors::DetID::ITS> + ; +// #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::trk::DPLDigitizerParam < o2::detectors::DetID::TRK>> + ; +// #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::trk::DPLDigitizerParam < o2::detectors::DetID::FT3>> + ; #endif diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/CMakeLists.txt b/Detectors/Upgrades/ALICE3/TRK/workflow/CMakeLists.txt new file mode 100644 index 0000000000000..c9f4099017717 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/CMakeLists.txt @@ -0,0 +1,32 @@ +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization +# or submit itself to any jurisdiction. + +o2_add_library(TRKWorkflow + SOURCES src/DigitReaderSpec.cxx + src/DigitWriterSpec.cxx + # src/RecoWorkflow.cxx + # src/ClusterWriterWorkflow.cxx + # src/ClustererSpec.cxx + # src/ClusterWriterSpec.cxx + # src/TrackerSpec.cxx + # src/TrackWriterSpec.cxx + # src/TrackReaderSpec.cxx + # src/VertexReaderSpec.cxx + PUBLIC_LINK_LIBRARIES O2::Framework + O2::SimConfig + O2::DataFormatsITSMFT + O2::SimulationDataFormat + O2::DPLUtils) + +# o2_add_executable(reco-workflow +# SOURCES src/trk-reco-workflow.cxx +# COMPONENT_NAME alice3-trk +# PUBLIC_LINK_LIBRARIES O2::TRKWorkflow) \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/DigitReaderSpec.h b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/DigitReaderSpec.h new file mode 100644 index 0000000000000..2a0acd792f4a9 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/DigitReaderSpec.h @@ -0,0 +1,87 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_TRK_DIGITREADER +#define O2_TRK_DIGITREADER + +#include "TFile.h" +#include "TTree.h" +#include "DataFormatsITSMFT/Digit.h" +#include "DataFormatsITSMFT/GBTCalibData.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" +#include "Headers/DataHeader.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "DetectorsCommonDataFormats/DetID.h" + +using namespace o2::framework; + +namespace o2 +{ +namespace trk +{ + +class DigitReader : public Task +{ + public: + DigitReader() = delete; + DigitReader(o2::detectors::DetID id, bool useMC, bool useCalib); + ~DigitReader() override = default; + void init(InitContext& ic) final; + void run(ProcessingContext& pc) final; + + protected: + void connectTree(const std::string& filename); + + std::vector mDigits, *mDigitsPtr = &mDigits; + std::vector mCalib, *mCalibPtr = &mCalib; + std::vector mDigROFRec, *mDigROFRecPtr = &mDigROFRec; + std::vector mDigMC2ROFs, *mDigMC2ROFsPtr = &mDigMC2ROFs; + + o2::header::DataOrigin mOrigin = o2::header::gDataOriginInvalid; + + std::unique_ptr mFile; + std::unique_ptr mTree; + + bool mUseMC = true; // use MC truth + bool mUseCalib = true; // send calib data + + std::string mDetName = ""; + std::string mDetNameLC = ""; + std::string mFileName = ""; + std::string mDigTreeName = "o2sim"; + std::string mDigitBranchName = "Digit"; + std::string mDigROFBranchName = "DigitROF"; + std::string mCalibBranchName = "Calib"; + + std::string mDigtMCTruthBranchName = "DigitMCTruth"; + std::string mDigtMC2ROFBranchName = "DigitMC2ROF"; +}; + +class TRKDigitReader : public DigitReader +{ + public: + TRKDigitReader(bool useMC = true, bool useCalib = false) + : DigitReader(o2::detectors::DetID::TRK, useMC, useCalib) + { + mOrigin = o2::header::gDataOriginTRK; + } +}; + +/// create a processor spec +/// read ITS/MFT Digit data from a root file +framework::DataProcessorSpec getTRKDigitReaderSpec(bool useMC = true, bool useCalib = false, std::string defname = "trkdigits.root"); + +} // namespace trk +} // namespace o2 + +#endif /* O2_TRK_DigitREADER */ diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/DigitWriterSpec.h b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/DigitWriterSpec.h new file mode 100644 index 0000000000000..9c37d4318bb0f --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/DigitWriterSpec.h @@ -0,0 +1,26 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef STEER_TRKDIGITWRITER_H_ +#define STEER_TRKDIGITWRITER_H_ + +#include "Framework/DataProcessorSpec.h" + +namespace o2 +{ +namespace trk +{ + +o2::framework::DataProcessorSpec getTRKDigitWriterSpec(bool mctruth = true, bool dec = false, bool calib = false); +} // namespace trk +} // end namespace o2 + +#endif /* STEER_TRKDIGITWRITER_H_ */ diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/src/DigitReaderSpec.cxx b/Detectors/Upgrades/ALICE3/TRK/workflow/src/DigitReaderSpec.cxx new file mode 100644 index 0000000000000..09bb1f12a48e4 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/src/DigitReaderSpec.cxx @@ -0,0 +1,139 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include + +#include "TTree.h" + +#include "Framework/ControlService.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/Logger.h" +#include "TRKWorkflow/DigitReaderSpec.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/ConstMCTruthContainer.h" +#include "SimulationDataFormat/IOMCTruthContainerView.h" +#include + +using namespace o2::framework; +using namespace o2::itsmft; + +namespace o2 +{ +namespace trk +{ + +DigitReader::DigitReader(o2::detectors::DetID id, bool useMC, bool useCalib) +{ + assert(id == o2::detectors::DetID::TRK); + mDetNameLC = mDetName = id.getName(); + mDigTreeName = "o2sim"; + + mDigitBranchName = mDetName + mDigitBranchName; + mDigROFBranchName = mDetName + mDigROFBranchName; + mCalibBranchName = mDetName + mCalibBranchName; + + mDigtMCTruthBranchName = mDetName + mDigtMCTruthBranchName; + mDigtMC2ROFBranchName = mDetName + mDigtMC2ROFBranchName; + + mUseMC = useMC; + mUseCalib = useCalib; + std::transform(mDetNameLC.begin(), mDetNameLC.end(), mDetNameLC.begin(), ::tolower); +} + +void DigitReader::init(InitContext& ic) +{ + mFileName = ic.options().get((mDetNameLC + "-digit-infile").c_str()); + connectTree(mFileName); +} + +void DigitReader::run(ProcessingContext& pc) +{ + auto ent = mTree->GetReadEntry() + 1; + assert(ent < mTree->GetEntries()); // this should not happen + + o2::dataformats::IOMCTruthContainerView* plabels = nullptr; + if (mUseMC) { + mTree->SetBranchAddress(mDigtMCTruthBranchName.c_str(), &plabels); + } + mTree->GetEntry(ent); + LOG(info) << mDetName << "DigitReader pushes " << mDigROFRec.size() << " ROFRecords, " + << mDigits.size() << " digits at entry " << ent; + + // This is a very ugly way of providing DataDescription, which anyway does not need to contain detector name. + // To be fixed once the names-definition class is ready + pc.outputs().snapshot(Output{mOrigin, "DIGITSROF", 0}, mDigROFRec); + pc.outputs().snapshot(Output{mOrigin, "DIGITS", 0}, mDigits); + if (mUseCalib) { + pc.outputs().snapshot(Output{mOrigin, "GBTCALIB", 0}, mCalib); + } + + if (mUseMC) { + auto& sharedlabels = pc.outputs().make>(Output{mOrigin, "DIGITSMCTR", 0}); + plabels->copyandflatten(sharedlabels); + delete plabels; + pc.outputs().snapshot(Output{mOrigin, "DIGITSMC2ROF", 0}, mDigMC2ROFs); + } + + if (mTree->GetReadEntry() + 1 >= mTree->GetEntries()) { + pc.services().get().endOfStream(); + pc.services().get().readyToQuit(QuitRequest::Me); + } +} + +void DigitReader::connectTree(const std::string& filename) +{ + mTree.reset(nullptr); // in case it was already loaded + mFile.reset(TFile::Open(filename.c_str())); + assert(mFile && !mFile->IsZombie()); + mTree.reset((TTree*)mFile->Get(mDigTreeName.c_str())); + assert(mTree); + + mTree->SetBranchAddress(mDigROFBranchName.c_str(), &mDigROFRecPtr); + mTree->SetBranchAddress(mDigitBranchName.c_str(), &mDigitsPtr); + if (mUseCalib) { + if (!mTree->GetBranch(mCalibBranchName.c_str())) { + throw std::runtime_error("GBT calibration data requested but not found in the tree"); + } + mTree->SetBranchAddress(mCalibBranchName.c_str(), &mCalibPtr); + } + if (mUseMC) { + if (!mTree->GetBranch(mDigtMC2ROFBranchName.c_str()) || !mTree->GetBranch(mDigtMCTruthBranchName.c_str())) { + throw std::runtime_error("MC data requested but not found in the tree"); + } + mTree->SetBranchAddress(mDigtMC2ROFBranchName.c_str(), &mDigMC2ROFsPtr); + } + LOG(info) << "Loaded tree from " << filename << " with " << mTree->GetEntries() << " entries"; +} + +DataProcessorSpec getTRKDigitReaderSpec(bool useMC, bool useCalib, std::string defname) +{ + std::vector outputSpec; + outputSpec.emplace_back("TRK", "DIGITS", 0, Lifetime::Timeframe); + outputSpec.emplace_back("TRK", "DIGITSROF", 0, Lifetime::Timeframe); + if (useCalib) { + outputSpec.emplace_back("TRK", "GBTCALIB", 0, Lifetime::Timeframe); + } + if (useMC) { + outputSpec.emplace_back("TRK", "DIGITSMCTR", 0, Lifetime::Timeframe); + outputSpec.emplace_back("TRK", "DIGITSMC2ROF", 0, Lifetime::Timeframe); + } + + return DataProcessorSpec{ + "trk-digit-reader", + Inputs{}, + outputSpec, + AlgorithmSpec{adaptFromTask(useMC, useCalib)}, + Options{ + {"trk-digit-infile", VariantType::String, defname, {"Name of the input digit file"}}}}; +} + +} // namespace trk +} // namespace o2 diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/src/DigitWriterSpec.cxx b/Detectors/Upgrades/ALICE3/TRK/workflow/src/DigitWriterSpec.cxx new file mode 100644 index 0000000000000..2a743551adddb --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/src/DigitWriterSpec.cxx @@ -0,0 +1,110 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @brief Processor spec for a ROOT file writer for ITSMFT digits + +#include "TRKWorkflow/DigitWriterSpec.h" +#include "DPLUtils/MakeRootTreeWriterSpec.h" +#include "DataFormatsITSMFT/Digit.h" +#include "DataFormatsITSMFT/GBTCalibData.h" +#include "Headers/DataHeader.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "SimulationDataFormat/ConstMCTruthContainer.h" +#include "SimulationDataFormat/IOMCTruthContainerView.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include +#include +#include + +using namespace o2::framework; +using SubSpecificationType = o2::framework::DataAllocator::SubSpecificationType; + +namespace o2 +{ +namespace trk +{ + +template +using BranchDefinition = MakeRootTreeWriterSpec::BranchDefinition; +using MCCont = o2::dataformats::ConstMCTruthContainer; + +/// create the processor spec +/// describing a processor receiving digits for ITS/MFT and writing them to file +DataProcessorSpec getDigitWriterSpec(bool mctruth, bool dec, bool calib, o2::header::DataOrigin detOrig, o2::detectors::DetID detId) +{ + std::string detStr = o2::detectors::DetID::getName(detId); + std::string detStrL = dec ? "o2_" : ""; // for decoded digits prepend by o2 + detStrL += detStr; + std::transform(detStrL.begin(), detStrL.end(), detStrL.begin(), ::tolower); + auto logger = [](std::vector const& inDigits) { + LOG(info) << "RECEIVED DIGITS SIZE " << inDigits.size(); + }; + + // the callback to be set as hook for custom action when the writer is closed + auto finishWriting = [](TFile* outputfile, TTree* outputtree) { + const auto* brArr = outputtree->GetListOfBranches(); + int64_t nent = 0; + for (const auto* brc : *brArr) { + int64_t n = ((const TBranch*)brc)->GetEntries(); + if (nent && (nent != n)) { + LOG(error) << "Branches have different number of entries"; + } + nent = n; + } + outputtree->SetEntries(nent); + outputtree->Write("", TObject::kOverwrite); + outputfile->Close(); + }; + + // handler for labels + // This is necessary since we can't store the original label buffer in a ROOT entry -- as is -- if it exceeds a certain size. + // We therefore convert it to a special split class. + auto fillLabels = [](TBranch& branch, std::vector const& labelbuffer, DataRef const& /*ref*/) { + o2::dataformats::ConstMCTruthContainerView labels(labelbuffer); + LOG(info) << "WRITING " << labels.getNElements() << " LABELS "; + + o2::dataformats::IOMCTruthContainerView outputcontainer; + auto ptr = &outputcontainer; + auto br = framework::RootTreeWriter::remapBranch(branch, &ptr); + outputcontainer.adopt(labelbuffer); + br->Fill(); + br->ResetAddress(); + }; + + return MakeRootTreeWriterSpec((detStr + "DigitWriter" + (dec ? "_dec" : "")).c_str(), + (detStrL + "digits.root").c_str(), + MakeRootTreeWriterSpec::TreeAttributes{"o2sim", "Digits tree"}, + MakeRootTreeWriterSpec::CustomClose(finishWriting), + // in case of labels we first read them as std::vector and process them correctly in the fillLabels hook + BranchDefinition>{InputSpec{"digitsMCTR", detOrig, "DIGITSMCTR", 0}, + (detStr + "DigitMCTruth").c_str(), + (mctruth ? 1 : 0), fillLabels}, + BranchDefinition>{InputSpec{"digitsMC2ROF", detOrig, "DIGITSMC2ROF", 0}, + (detStr + "DigitMC2ROF").c_str(), + (mctruth ? 1 : 0)}, + BranchDefinition>{InputSpec{"digits", detOrig, "DIGITS", 0}, + (detStr + "Digit").c_str(), + logger}, + BranchDefinition>{InputSpec{"calib", detOrig, "GBTCALIB", 0}, + (detStr + "Calib").c_str(), + (calib ? 1 : 0)}, + BranchDefinition>{InputSpec{"digitsROF", detOrig, "DIGITSROF", 0}, + (detStr + "DigitROF").c_str()})(); +} + +DataProcessorSpec getTRKDigitWriterSpec(bool mctruth, bool dec, bool calib) +{ + return getDigitWriterSpec(mctruth, dec, calib, o2::header::gDataOriginTRK, o2::detectors::DetID::TRK); +} + +} // end namespace trk +} // end namespace o2 diff --git a/Steer/DigitizerWorkflow/CMakeLists.txt b/Steer/DigitizerWorkflow/CMakeLists.txt index 1b839ba462b63..babc5fce4d864 100644 --- a/Steer/DigitizerWorkflow/CMakeLists.txt +++ b/Steer/DigitizerWorkflow/CMakeLists.txt @@ -29,6 +29,7 @@ o2_add_executable(digitizer-workflow src/ZDCDigitizerSpec.cxx src/TOFDigitizerSpec.cxx $<$:src/ITS3DigitizerSpec.cxx> + $<$:src/TRKDigitizerSpec.cxx> PUBLIC_LINK_LIBRARIES O2::Framework O2::Steer O2::CommonConstants @@ -67,7 +68,9 @@ o2_add_executable(digitizer-workflow O2::DetectorsRaw $<$:O2::ITS3Simulation> $<$:O2::ITS3Workflow> - $<$:O2::ITS3Align>) + $<$:O2::ITS3Align> + $<$:O2::TRKSimulation> + $<$:O2::TRKWorkflow>) o2_add_executable(mctruth-testworkflow diff --git a/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx b/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx index a30294a240fb0..75141425f7c49 100644 --- a/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx +++ b/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx @@ -44,6 +44,10 @@ // for ITS3 #include "ITS3DigitizerSpec.h" #include "ITS3Workflow/DigitWriterSpec.h" + +// for alice 3 TRK +#include "TRKDigitizerSpec.h" +#include "TRKWorkflow/DigitWriterSpec.h" #endif // for TOF @@ -632,6 +636,15 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) // // connect ITS digit writer specs.emplace_back(o2::its3::getITS3DigitWriterSpec(mctruth)); } + + // the ALICE 3 TRK part + if (isEnabled(o2::detectors::DetID::TRK)) { + detList.emplace_back(o2::detectors::DetID::TRK); + // connect the ALICE 3 TRK digitization + specs.emplace_back(o2::trk::getTRKDigitizerSpec(fanoutsize++, mctruth)); + // connect the ALICE 3 TRK digit writer + specs.emplace_back(o2::trk::getTRKDigitWriterSpec(mctruth)); + } #endif // the MFT part diff --git a/Steer/DigitizerWorkflow/src/TRKDigitizerSpec.cxx b/Steer/DigitizerWorkflow/src/TRKDigitizerSpec.cxx new file mode 100644 index 0000000000000..1d923e6c3d4d2 --- /dev/null +++ b/Steer/DigitizerWorkflow/src/TRKDigitizerSpec.cxx @@ -0,0 +1,304 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "TRKDigitizerSpec.h" +#include "Framework/ControlService.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/CCDBParamSpec.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/DataRefUtils.h" +#include "Framework/Lifetime.h" +#include "Framework/Task.h" +#include "Steer/HitProcessingManager.h" +#include "DataFormatsITSMFT/Digit.h" +#include "SimulationDataFormat/ConstMCTruthContainer.h" +#include "DetectorsBase/BaseDPLDigitizer.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "DetectorsCommonDataFormats/SimTraits.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "TRKSimulation/Digitizer.h" +#include "TRKSimulation/DPLDigitizerParam.h" +#include "ITSMFTBase/DPLAlpideParam.h" +#include "TRKBase/GeometryTGeo.h" +#include "TRKBase/TRKBaseParam.h" +// #includTRKAlign/MisalignmentManager.h" + +#include +#include + +#include + +using namespace o2::framework; +using SubSpecificationType = o2::framework::DataAllocator::SubSpecificationType; + +namespace +{ +std::vector makeOutChannels(o2::header::DataOrigin detOrig, bool mctruth) +{ + std::vector outputs; + outputs.emplace_back(detOrig, "DIGITS", 0, Lifetime::Timeframe); + outputs.emplace_back(detOrig, "DIGITSROF", 0, Lifetime::Timeframe); + if (mctruth) { + outputs.emplace_back(detOrig, "DIGITSMC2ROF", 0, Lifetime::Timeframe); + outputs.emplace_back(detOrig, "DIGITSMCTR", 0, Lifetime::Timeframe); + } + outputs.emplace_back(detOrig, "ROMode", 0, Lifetime::Timeframe); + return outputs; +} +} // namespace + +namespace o2::trk +{ +using namespace o2::base; +class TRKDPLDigitizerTask : BaseDPLDigitizer +{ + public: + using BaseDPLDigitizer::init; + + TRKDPLDigitizerTask(bool mctruth = true) : BaseDPLDigitizer(InitServices::FIELD | InitServices::GEOM), mWithMCTruth(mctruth) {} + + void initDigitizerTask(framework::InitContext& ic) override + { + mDisableQED = ic.options().get("disable-qed"); + } + + void run(framework::ProcessingContext& pc) + { + if (mFinished) { + return; + } + updateTimeDependentParams(pc); + + // read collision context from input + auto context = pc.inputs().get("collisioncontext"); + context->initSimChains(mID, mSimChains); + const bool withQED = context->isQEDProvided() && !mDisableQED; + auto& timesview = context->getEventRecords(withQED); + LOG(info) << "GOT " << timesview.size() << " COLLISION TIMES"; + LOG(info) << "SIMCHAINS " << mSimChains.size(); + + // if there is nothing to do ... return + if (timesview.empty()) { + return; + } + TStopwatch timer; + timer.Start(); + LOG(info) << " CALLING TRK DIGITIZATION "; + + // mDigitizer.setDigits(&mDigits); + mDigitizer.setROFRecords(&mROFRecords); + mDigitizer.setMCLabels(&mLabels); + + // digits are directly put into DPL owned resource + auto& digitsAccum = pc.outputs().make>(Output{mOrigin, "DIGITS", 0}); + + auto accumulate = [this, &digitsAccum]() { + // accumulate result of single event processing, called after processing every event supplied + // AND after the final flushing via digitizer::fillOutputContainer + if (mDigits.empty()) { + return; // no digits were flushed, nothing to accumulate + } + auto ndigAcc = digitsAccum.size(); + std::copy(mDigits.begin(), mDigits.end(), std::back_inserter(digitsAccum)); + + // fix ROFrecords references on ROF entries + auto nROFRecsOld = mROFRecordsAccum.size(); + + for (int i = 0; i < mROFRecords.size(); i++) { + auto& rof = mROFRecords[i]; + rof.setFirstEntry(ndigAcc + rof.getFirstEntry()); + rof.print(); + + if (mFixMC2ROF < mMC2ROFRecordsAccum.size()) { // fix ROFRecord entry in MC2ROF records + for (int m2rid = mFixMC2ROF; m2rid < mMC2ROFRecordsAccum.size(); m2rid++) { + // need to register the ROFRecors entry for MC event starting from this entry + auto& mc2rof = mMC2ROFRecordsAccum[m2rid]; + if (rof.getROFrame() == mc2rof.minROF) { + mFixMC2ROF++; + mc2rof.rofRecordID = nROFRecsOld + i; + mc2rof.print(); + } + } + } + } + + std::copy(mROFRecords.begin(), mROFRecords.end(), std::back_inserter(mROFRecordsAccum)); + if (mWithMCTruth) { + mLabelsAccum.mergeAtBack(mLabels); + } + LOG(info) << "Added " << mDigits.size() << " digits "; + // clean containers from already accumulated stuff + mLabels.clear(); + mDigits.clear(); + mROFRecords.clear(); + }; // and accumulate lambda + + auto& eventParts = context->getEventParts(withQED); + // loop over all composite collisions given from context (aka loop over all the interaction records) + const int bcShift = mDigitizer.getParams().getROFrameBiasInBC(); + // loop over all composite collisions given from context (aka loop over all the interaction records) + for (size_t collID = 0; collID < timesview.size(); ++collID) { + auto irt = timesview[collID]; + if (irt.toLong() < bcShift) { // due to the ROF misalignment the collision would go to negative ROF ID, discard + continue; + } + irt -= bcShift; // account for the ROF start shift + + mDigitizer.setEventTime(irt); + mDigitizer.resetEventROFrames(); // to estimate min/max ROF for this collID + // for each collision, loop over the constituents event and source IDs + // (background signal merging is basically taking place here) + for (auto& part : eventParts[collID]) { + + // get the hits for this event and this source + mHits.clear(); + context->retrieveHits(mSimChains, o2::detectors::SimTraits::DETECTORBRANCHNAMES[mID][0].c_str(), part.sourceID, part.entryID, &mHits); + + if (!mHits.empty()) { + LOG(debug) << "For collision " << collID << " eventID " << part.entryID + << " found " << mHits.size() << " hits "; + mDigitizer.process(&mHits, part.entryID, part.sourceID); // call actual digitization procedure + } + } + mMC2ROFRecordsAccum.emplace_back(collID, -1, mDigitizer.getEventROFrameMin(), mDigitizer.getEventROFrameMax()); + accumulate(); + } + mDigitizer.fillOutputContainer(); + accumulate(); + + // here we have all digits and labels and we can send them to consumer (aka snapshot it onto output) + + pc.outputs().snapshot(Output{mOrigin, "DIGITSROF", 0}, mROFRecordsAccum); + if (mWithMCTruth) { + pc.outputs().snapshot(Output{mOrigin, "DIGITSMC2ROF", 0}, mMC2ROFRecordsAccum); + auto& sharedlabels = pc.outputs().make>(Output{mOrigin, "DIGITSMCTR", 0}); + mLabelsAccum.flatten_to(sharedlabels); + // free space of existing label containers + mLabels.clear_andfreememory(); + mLabelsAccum.clear_andfreememory(); + } + LOG(info) << mID.getName() << ": Sending ROMode= " << mROMode << " to GRPUpdater"; + pc.outputs().snapshot(Output{mOrigin, "ROMode", 0}, mROMode); + + timer.Stop(); + LOG(info) << "Digitization took " << timer.CpuTime() << "s"; + + // we should be only called once; tell DPL that this process is ready to exit + pc.services().get().readyToQuit(QuitRequest::Me); + + mFinished = true; + } + + void updateTimeDependentParams(ProcessingContext& pc) + { + static bool initOnce{false}; + if (!initOnce) { + initOnce = true; + auto& digipar = mDigitizer.getParams(); + + // configure digitizer + o2::trk::GeometryTGeo* geom = o2::trk::GeometryTGeo::Instance(); + geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); // make sure L2G matrices are loaded + mDigitizer.setGeometry(geom); + + const auto& dopt = o2::trk::DPLDigitizerParam::Instance(); + pc.inputs().get*>("ITS_alppar"); + const auto& aopt = o2::itsmft::DPLAlpideParam::Instance(); + digipar.setContinuous(dopt.continuous); + digipar.setROFrameBiasInBC(aopt.roFrameBiasInBC); + if (dopt.continuous) { + auto frameNS = aopt.roFrameLengthInBC * o2::constants::lhc::LHCBunchSpacingNS; + digipar.setROFrameLengthInBC(aopt.roFrameLengthInBC); + digipar.setROFrameLength(frameNS); // RO frame in ns + digipar.setStrobeDelay(aopt.strobeDelay); // Strobe delay wrt beginning of the RO frame, in ns + digipar.setStrobeLength(aopt.strobeLengthCont > 0 ? aopt.strobeLengthCont : frameNS - aopt.strobeDelay); // Strobe length in ns + } else { + digipar.setROFrameLength(aopt.roFrameLengthTrig); // RO frame in ns + digipar.setStrobeDelay(aopt.strobeDelay); // Strobe delay wrt beginning of the RO frame, in ns + digipar.setStrobeLength(aopt.strobeLengthTrig); // Strobe length in ns + } + // parameters of signal time response: flat-top duration, max rise time and q @ which rise time is 0 + digipar.getSignalShape().setParameters(dopt.strobeFlatTop, dopt.strobeMaxRiseTime, dopt.strobeQRiseTime0); + digipar.setChargeThreshold(dopt.chargeThreshold); // charge threshold in electrons + digipar.setNoisePerPixel(dopt.noisePerPixel); // noise level + digipar.setTimeOffset(dopt.timeOffset); + digipar.setNSimSteps(dopt.nSimSteps); + + mROMode = digipar.isContinuous() ? o2::parameters::GRPObject::CONTINUOUS : o2::parameters::GRPObject::PRESENT; + LOG(info) << mID.getName() << " simulated in " + << ((mROMode == o2::parameters::GRPObject::CONTINUOUS) ? "CONTINUOUS" : "TRIGGERED") + << " RO mode"; + + // if (oTRKParams::Instance().useDeadChannelMap) { + // pc.inputs().get("TRK_dead"); // trigger final ccdb update + // } + + // init digitizer + mDigitizer.init(); + } + // Other time-dependent parameters can be added below + } + + void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) + { + if (matcher == ConcreteDataMatcher(detectors::DetID::ITS, "ALPIDEPARAM", 0)) { + LOG(info) << mID.getName() << " Alpide param updated"; + const auto& par = o2::itsmft::DPLAlpideParam::Instance(); + par.printKeyValues(); + return; + } + // if (matcher == ConcreteDataMatcher(mOrigin, "DEADMAP", 0)) { + // LOG(info) << mID.getName() << " static dead map updated"; + // mDigitizer.setDeadChannelsMap((o2::itsmft::NoiseMap*)obj); + // return; + // } + } + + private: + bool mWithMCTruth{true}; + bool mFinished{false}; + bool mDisableQED{false}; + const o2::detectors::DetID mID{o2::detectors::DetID::TRK}; + const o2::header::DataOrigin mOrigin{o2::header::gDataOriginTRK}; + o2::trk::Digitizer mDigitizer{}; + std::vector mDigits{}; + std::vector mROFRecords{}; + std::vector mROFRecordsAccum{}; + std::vector mHits{}; + std::vector* mHitsP{&mHits}; + o2::dataformats::MCTruthContainer mLabels{}; + o2::dataformats::MCTruthContainer mLabelsAccum{}; + std::vector mMC2ROFRecordsAccum{}; + std::vector mSimChains{}; + + int mFixMC2ROF = 0; // 1st entry in mc2rofRecordsAccum to be fixed for ROFRecordID + o2::parameters::GRPObject::ROMode mROMode = o2::parameters::GRPObject::PRESENT; // readout mode +}; + +DataProcessorSpec getTRKDigitizerSpec(int channel, bool mctruth) +{ + std::string detStr = o2::detectors::DetID::getName(o2::detectors::DetID::TRK); + auto detOrig = o2::header::gDataOriginTRK; + std::vector inputs; + inputs.emplace_back("collisioncontext", "SIM", "COLLISIONCONTEXT", static_cast(channel), Lifetime::Timeframe); + inputs.emplace_back("ITS_alppar", "ITS", "ALPIDEPARAM", 0, Lifetime::Condition, ccdbParamSpec("ITS/Config/AlpideParam")); + // if (oTRKParams::Instance().useDeadChannelMap) { + // inputs.emplace_back("TRK_dead", "TRK", "DEADMAP", 0, Lifetime::Condition, ccdbParamSpec("TRK/Calib/DeadMap")); + // } + + return DataProcessorSpec{detStr + "Digitizer", + inputs, makeOutChannels(detOrig, mctruth), + AlgorithmSpec{adaptFromTask(mctruth)}, + Options{{"disable-qed", o2::framework::VariantType::Bool, false, {"disable QED handling"}}}}; +} + +} // namespace o2::trk diff --git a/Steer/DigitizerWorkflow/src/TRKDigitizerSpec.h b/Steer/DigitizerWorkflow/src/TRKDigitizerSpec.h new file mode 100644 index 0000000000000..5a1a59c3b9f5e --- /dev/null +++ b/Steer/DigitizerWorkflow/src/TRKDigitizerSpec.h @@ -0,0 +1,24 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef STEER_DIGITIZERWORKFLOW_TRKDIGITIZER_H_ +#define STEER_DIGITIZERWORKFLOW_TRKDIGITIZER_H_ + +#include "Framework/DataProcessorSpec.h" + +namespace o2::trk +{ +o2::framework::DataProcessorSpec getTRKDigitizerSpec(int channel, bool mctruth = true); +} +// namespace o2::trk +// end namespace o2 + +#endif