diff --git a/Detectors/HMPID/base/include/HMPIDBase/Digit.h b/Detectors/HMPID/base/include/HMPIDBase/Digit.h index a19aeb33b2f7e..381f2556b273b 100644 --- a/Detectors/HMPID/base/include/HMPIDBase/Digit.h +++ b/Detectors/HMPID/base/include/HMPIDBase/Digit.h @@ -31,48 +31,60 @@ class Digit : public DigitBase Digit(float charge) : mQ(charge) {} float getCharge() const { return mQ; } int getPadID() const { return mPad; } + // convenience conversion to x-y pad coordinates + int getPx() const { return Param::A2X(mPad); } + int getPy() const { return Param::A2Y(mPad); } - Digit(HitType const& hit) + static void getPadAndTotalCharge(HitType const& hit, int& chamber, int& pc, int& px, int& py, float& totalcharge) { - int pc, px, py; float localX; float localY; - - // chamber number is in detID - const auto chamber = hit.GetDetectorID(); - + chamber = hit.GetDetectorID(); double tmp[3] = { hit.GetX(), hit.GetY(), hit.GetZ() }; Param::Instance()->Mars2Lors(chamber, tmp, localX, localY); - Param::Lors2Pad(localX, localY, pc, px, py); - // TODO: check if this digit is valid - // mark as invalid otherwise - - // calculate pad id - mPad = Param::Abs(chamber, pc, px, py); + totalcharge = Digit::QdcTot(hit.GetEnergyLoss(), hit.GetTime(), pc, px, py, localX, localY); + } - // calculate charge - mQ = /*fQHit **/ InMathieson(localX, localY); + static float getFractionalContributionForPad(HitType const& hit, int somepad) + { + float localX; + float localY; - // what about time stamp?? + // chamber number is in detID + const auto chamber = hit.GetDetectorID(); + double tmp[3] = { hit.GetX(), hit.GetY(), hit.GetZ() }; + // converting chamber id and hit coordiates to local coordinates + Param::Instance()->Mars2Lors(chamber, tmp, localX, localY); + // calculate charge fraction in given pad + return Digit::InMathieson(localX, localY, somepad); } + Digit(int pad, float charge) : mPad(pad), mQ(charge) {} + + // add charge to existing digit + void addCharge(float q) { mQ += q; } + private: float mQ = 0.; int mPad = 0.; // -1 indicates invalid digit - float LorsX() const { return Param::LorsX(Param::A2P(mPad), Param::A2X(mPad)); } //center of the pad x, [cm] - float LorsY() const { return Param::LorsY(Param::A2P(mPad), Param::A2Y(mPad)); } //center of the pad y, [cm] + static float LorsX(int pad) { return Param::LorsX(Param::A2P(pad), Param::A2X(pad)); } //center of the pad x, [cm] + static float LorsY(int pad) { return Param::LorsY(Param::A2P(pad), Param::A2Y(pad)); } //center of the pad y, [cm] + + // determines the total charge created by a hit + // might modify the localX, localY coordiates associated to the hit + static float QdcTot(float e, float time, int pc, int px, int py, float& localX, float& localY); - float IntPartMathiX(float x) const + static float IntPartMathiX(float x, int pad) { // Integration of Mathieson. // This is the answer to electrostatic problem of charge distrubution in MWPC described elsewhere. (NIM A370(1988)602-603) // Arguments: x,y- position of the center of Mathieson distribution // Returns: a charge fraction [0-1] imposed into the pad - auto shift1 = -LorsX() + 0.5 * Param::SizePadX(); - auto shift2 = -LorsX() - 0.5 * Param::SizePadX(); + auto shift1 = -LorsX(pad) + 0.5 * Param::SizePadX(); + auto shift2 = -LorsX(pad) - 0.5 * Param::SizePadX(); auto ux1 = Param::SqrtK3x() * TMath::TanH(Param::K2x() * (x + shift1) / Param::PitchAnodeCathode()); auto ux2 = Param::SqrtK3x() * TMath::TanH(Param::K2x() * (x + shift2) / o2::hmpid::Param::PitchAnodeCathode()); @@ -81,14 +93,14 @@ class Digit : public DigitBase } //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - Double_t IntPartMathiY(Double_t y) const + static Double_t IntPartMathiY(Double_t y, int pad) { // Integration of Mathieson. // This is the answer to electrostatic problem of charge distrubution in MWPC described elsewhere. (NIM A370(1988)602-603) // Arguments: x,y- position of the center of Mathieson distribution // Returns: a charge fraction [0-1] imposed into the pad - Double_t shift1 = -LorsY() + 0.5 * o2::hmpid::Param::SizePadY(); - Double_t shift2 = -LorsY() - 0.5 * o2::hmpid::Param::SizePadY(); + Double_t shift1 = -LorsY(pad) + 0.5 * o2::hmpid::Param::SizePadY(); + Double_t shift2 = -LorsY(pad) - 0.5 * o2::hmpid::Param::SizePadY(); Double_t uy1 = Param::SqrtK3y() * TMath::TanH(Param::K2y() * (y + shift1) / Param::PitchAnodeCathode()); Double_t uy2 = Param::SqrtK3y() * TMath::TanH(Param::K2y() * (y + shift2) / Param::PitchAnodeCathode()); @@ -96,12 +108,12 @@ class Digit : public DigitBase return Param::K4y() * (TMath::ATan(uy2) - TMath::ATan(uy1)); } - float InMathieson(float localX, float localY) const + static float InMathieson(float localX, float localY, int pad) { - return 4. * IntPartMathiX(localX) * IntPartMathiY(localY); + return 4. * Digit::IntPartMathiX(localX, pad) * Digit::IntPartMathiY(localY, pad); } - ClassDefNV(Digit, 1); + ClassDefNV(Digit, 2); }; } // namespace hmpid diff --git a/Detectors/HMPID/base/src/Digit.cxx b/Detectors/HMPID/base/src/Digit.cxx index 6f21fa4c4d623..de4c9cea63ff4 100644 --- a/Detectors/HMPID/base/src/Digit.cxx +++ b/Detectors/HMPID/base/src/Digit.cxx @@ -9,8 +9,43 @@ // or submit itself to any jurisdiction. #include "HMPIDBase/Digit.h" -#include +#include "HMPIDBase/Param.h" +#include "TRandom.h" +#include "TMath.h" using namespace o2::hmpid; ClassImp(o2::hmpid::Digit); + +float Digit::QdcTot(float e, float time, int pc, int px, int py, float& localX, float& localY) +{ + // Samples total charge associated to a hit + // Arguments: e- hit energy [GeV] for mip Eloss for photon Etot + // Returns: total QDC + float Q = 0; + if (time > 1.2e-6) { + Q = 0; + } + if (py < 0) { + return 0; + } else { + float y = Param::LorsY(pc, py); + localY = ((y - localY) > 0) ? y - 0.2 : y + 0.2; //shift to the nearest anod wire + + float x = (localX > 66.6) ? localX - 66.6 : localX; //sagita is for PC (0-64) and not for chamber + float qdcEle = 34.06311 + 0.2337070 * x + 5.807476e-3 * x * x - 2.956471e-04 * x * x * x + 2.310001e-06 * x * x * x * x; //reparametrised from DiMauro + + int iNele = int((e / 26e-9) * 0.8); + if (iNele < 1) { + iNele = 1; //number of electrons created by hit, if photon e=0 implies iNele=1 + } + for (Int_t i = 1; i <= iNele; i++) { + double rnd = gRandom->Rndm(); + if (rnd == 0) { + rnd = 1e-12; //1e-12 is a protection against 0 from rndm + } + Q -= qdcEle * TMath::Log(rnd); + } + } + return Q; +} diff --git a/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h b/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h index 1c3462fb6ddba..38b26e4870d40 100644 --- a/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h +++ b/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h @@ -13,6 +13,8 @@ #include "HMPIDBase/Digit.h" #include "HMPIDSimulation/Detector.h" // for the hit +#include "SimulationDataFormat/MCTruthContainer.h" +#include "SimulationDataFormat/MCCompLabel.h" #include namespace o2 @@ -27,10 +29,24 @@ class HMPIDDigitizer void setEventID(int eventID) { mEventID = eventID; } void setSrcID(int sID) { mSrcID = sID; } + // user can pass a label container to be filled ... this activates the label mechanism + // the passed label container can be readout after call to process + void setLabelContainer(o2::dataformats::MCTruthContainer* labels) + { + mRegisteredLabelContainer = labels; + } + // this will process hits and fill the digit vector with digits which are finalized void process(std::vector const&, std::vector& digit); private: + void zeroSuppress(std::vector const& digits, std::vector& newdigits, + o2::dataformats::MCTruthContainer const& labels, + o2::dataformats::MCTruthContainer* newlabels); + + float getThreshold(o2::hmpid::Digit const&) const; // gives back threshold to apply for a certain digit + // (using noise and other tables for pad) + double mTime = 0.; int mEventID = 0; int mSrcID = 0; @@ -39,7 +55,18 @@ class HMPIDDigitizer std::vector mSummable; std::vector mFinal; - // other stuff needed for digitizaton + std::vector mDigits; // internal store for digits + + //static constexpr int HMPID_NUMBEROFPADS = 161280; + //std::array mIndexForPad = { -1 }; //! mapping of pad to digit index + + std::map mIndexForPad; //! logarithmic mapping of pad to digit index + + std::vector mInvolvedPads; //! list of pads where digits created + + // other stuff needed for digitization + o2::dataformats::MCTruthContainer mTmpLabelContainer; // temp label container as workspace + o2::dataformats::MCTruthContainer* mRegisteredLabelContainer = nullptr; // label container to be filled ClassDefNV(HMPIDDigitizer, 1); }; diff --git a/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx b/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx index d4d785d724ec1..925efc4f1bbf9 100644 --- a/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx +++ b/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx @@ -15,13 +15,99 @@ using namespace o2::hmpid; ClassImp(HMPIDDigitizer); +float HMPIDDigitizer::getThreshold(o2::hmpid::Digit const& digiti) const +{ + // TODO: implement like in AliRoot some thresholding depending on conditions ... + return 4.; +} + +// applies threshold to digits; removes the ones below a certain charge threshold +void HMPIDDigitizer::zeroSuppress(std::vector const& digits, std::vector& newdigits, + o2::dataformats::MCTruthContainer const& labels, + o2::dataformats::MCTruthContainer* newlabels) +{ + int index = 0; + for (auto& digit : digits) { + if (digit.getCharge() >= getThreshold(digit)) { + newdigits.push_back(digit); + + if (newlabels) { + // copy the labels to the new place with the right new index + newlabels->addElements(newdigits.size() - 1, labels.getLabels(index)); + } + } + index++; + } +} + // this will process hits and fill the digit vector with digits which are finalized void HMPIDDigitizer::process(std::vector const& hits, std::vector& digits) { - // this is a very simple variant that creates one digit from one hit - // conversion is done in the digit constructor - // TODO: introduce cross-talk, pile-up, etc. + mIndexForPad.clear(); + mInvolvedPads.clear(); + for (auto& hit : hits) { - digits.emplace_back(hit); + int chamber, pc, px, py; + float totalQ; + // retrieves center pad and the total charge + Digit::getPadAndTotalCharge(hit, chamber, pc, px, py, totalQ); + + if (px < 0 || py < 0) { + continue; + } + + // determine which pads to loop over + std::array allpads; + int counter = 0; + for (int nx = -1; nx <= 1; ++nx) { + for (int ny = -1; ny <= 1; ++ny) { + allpads[counter] = Param::Abs(chamber, pc, px + nx, py + ny); + counter++; + } + } + + for (auto& pad : allpads) { + auto iter = mIndexForPad.find(pad); + int index = -1; + if (iter != mIndexForPad.end()) { + index = iter->second; + } + // auto index = mIndexForPad[pad]; + float fraction = Digit::getFractionalContributionForPad(hit, pad); + // LOG(INFO) << "FRACTION ON PAD " << pad << " IS " << fraction; + if (index != -1) { + // digit exists ... reuse + auto& digit = mDigits[index]; + digit.addCharge(totalQ * fraction); + + if (mRegisteredLabelContainer) { + auto labels = mTmpLabelContainer.getLabels(index); + o2::MCCompLabel newlabel(hit.GetTrackID(), mEventID, mSrcID); + bool newlabelneeded = true; + for (auto& l : labels) { + if (l == newlabel) { + newlabelneeded = false; + break; + } + } + if (newlabelneeded) { + mTmpLabelContainer.addElementRandomAccess(index, newlabel); + } + } + } else { + // create digit ... and register + mDigits.emplace_back(pad, totalQ * fraction); + mIndexForPad[pad] = mDigits.size() - 1; + mInvolvedPads.emplace_back(pad); + + if (mRegisteredLabelContainer) { + // add label for this digit + mTmpLabelContainer.addElement(mDigits.size() - 1, o2::MCCompLabel(hit.GetTrackID(), mEventID, mSrcID)); + } + } + } } + + // apply zero suppression (removes digits in the noise) + zeroSuppress(mDigits, digits, mTmpLabelContainer, mRegisteredLabelContainer); } diff --git a/Steer/DigitizerWorkflow/CMakeLists.txt b/Steer/DigitizerWorkflow/CMakeLists.txt index 03d41f1c1b43b..22a4217e5fd43 100644 --- a/Steer/DigitizerWorkflow/CMakeLists.txt +++ b/Steer/DigitizerWorkflow/CMakeLists.txt @@ -37,7 +37,6 @@ O2_GENERATE_EXECUTABLE( src/EMCALDigitizerSpec.cxx src/EMCALDigitWriterSpec.cxx src/HMPIDDigitizerSpec.cxx - src/HMPIDDigitWriterSpec.cxx src/MCHDigitizerSpec.cxx src/TRDDigitizerSpec.cxx src/GRPUpdaterSpec.cxx diff --git a/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.cxx b/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.cxx deleted file mode 100644 index bd62e7c44b61f..0000000000000 --- a/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.cxx +++ /dev/null @@ -1,73 +0,0 @@ -// Copyright CERN and copyright holders of ALICE O2. This software is -// distributed under the terms of the GNU General Public License v3 (GPL -// Version 3), copied verbatim in the file "COPYING". -// -// See http://alice-o2.web.cern.ch/license for full licensing information. -// -// 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 "HMPIDDigitWriterSpec.h" -#include "Framework/ControlService.h" -#include "Framework/DataProcessorSpec.h" -#include "Framework/DataRefUtils.h" -#include "Framework/Lifetime.h" -#include "Framework/Task.h" -#include "HMPIDBase/Digit.h" - -using namespace o2::framework; -using SubSpecificationType = o2::framework::DataAllocator::SubSpecificationType; - -namespace o2 -{ -namespace hmpid -{ - -class HMPIDDPLDigitWriterTask -{ - public: - void init(framework::InitContext& ic) - { - LOG(INFO) << "initializing HMPID digit writer"; - } - - void run(framework::ProcessingContext& pc) - { - static bool finished = false; - if (finished) { - return; - } - - LOG(INFO) << "Doing HMPID digit IO"; - auto digits = pc.inputs().get*>("hmpiddigits"); - LOG(INFO) << "HMPID received " << digits->size() << " digits"; - for (auto& digit : *digits) { - LOG(INFO) << "Have digit with charge " << digit.getCharge(); - } - - // we should be only called once; tell DPL that this process is ready to exit - pc.services().get().readyToQuit(false); - finished = true; - } -}; - -o2::framework::DataProcessorSpec getHMPIDDigitWriterSpec() -{ - // create the full data processor spec using - // a name identifier - // input description - // algorithmic description (here a lambda getting called once to setup the actual processing function) - // options that can be used for this processor (here: input file names where to take the hits) - return DataProcessorSpec{ - "HMPIDDigitWriternn", - Inputs{ InputSpec{ "hmpiddigits", "HMP", "DIGITS", 0, Lifetime::Timeframe } }, - {}, - AlgorithmSpec{ adaptFromTask() }, - Options{ - { "digitFile", VariantType::String, "hmpiddigits.root", { "filename for HMPID digits" } } } - }; -} - -} // end namespace hmpid -} // end namespace o2 diff --git a/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.h b/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.h index c03355d580160..4b26b7d60fe49 100644 --- a/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.h +++ b/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.h @@ -8,17 +8,35 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -#ifndef STEER_DIGITIZERWORKFLOW_SRC_HMPIDDIGITWRITERSPEC_H_ -#define STEER_DIGITIZERWORKFLOW_SRC_HMPIDDIGITWRITERSPEC_H_ +#ifndef STEER_DIGITIZERWORKFLOW_SRC_HMPDIGITWRITERSPEC_H_ +#define STEER_DIGITIZERWORKFLOW_SRC_HMPDIGITWRITERSPEC_H_ #include "Framework/DataProcessorSpec.h" +#include "Utils/MakeRootTreeWriterSpec.h" +#include "Framework/InputSpec.h" +#include "HMPIDBase/Digit.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "SimulationDataFormat/MCCompLabel.h" namespace o2 { namespace hmpid { -o2::framework::DataProcessorSpec getHMPIDDigitWriterSpec(); +template +using BranchDefinition = framework::MakeRootTreeWriterSpec::BranchDefinition; + +o2::framework::DataProcessorSpec getHMPIDDigitWriterSpec() +{ + using InputSpec = framework::InputSpec; + using MakeRootTreeWriterSpec = framework::MakeRootTreeWriterSpec; + return MakeRootTreeWriterSpec("HMPDigitWriter", + "hmpiddigits.root", + "o2sim", + 1, + BranchDefinition>{ InputSpec{ "digitinput", "HMP", "DIGITS" }, "HMPDigit" }, + BranchDefinition>{ InputSpec{ "labelinput", "HMP", "DIGITLBL" }, "HMPDigitLabels" })(); +} } // end namespace hmpid } // end namespace o2 diff --git a/Steer/DigitizerWorkflow/src/HMPIDDigitizerSpec.cxx b/Steer/DigitizerWorkflow/src/HMPIDDigitizerSpec.cxx index 5225fb0cc59c3..06e6771b10387 100644 --- a/Steer/DigitizerWorkflow/src/HMPIDDigitizerSpec.cxx +++ b/Steer/DigitizerWorkflow/src/HMPIDDigitizerSpec.cxx @@ -93,6 +93,7 @@ class HMPIDDPLDigitizerTask auto& eventParts = context->getEventParts(); std::vector digitsAccum; // accumulator for digits + o2::dataformats::MCTruthContainer labelAccum; // timeframe accumulator for labels // loop over all composite collisions given from context // (aka loop over all the interaction records) @@ -111,17 +112,18 @@ class HMPIDDPLDigitizerTask LOG(INFO) << "For collision " << collID << " eventID " << part.entryID << " found HMP " << hits.size() << " hits "; std::vector digits; // digits which get filled + o2::dataformats::MCTruthContainer labels; // labels which get filled + mDigitizer.setLabelContainer(&labels); mDigitizer.process(hits, digits); LOG(INFO) << "HMPID obtained " << digits.size() << " digits "; - for (auto& d : digits) { - LOG(INFO) << "CHARGE " << d.getCharge(); - LOG(INFO) << "PAD " << d.getPadID(); - } + LOG(INFO) << "NUMBER OF LABEL OBTAINED " << labels.getNElements(); std::copy(digits.begin(), digits.end(), std::back_inserter(digitsAccum)); + labelAccum.mergeAtBack(labels); } } pc.outputs().snapshot(Output{ "HMP", "DIGITS", 0, Lifetime::Timeframe }, digitsAccum); + pc.outputs().snapshot(Output{ "HMP", "DIGITLBL", 0, Lifetime::Timeframe }, labelAccum); LOG(INFO) << "HMP: Sending ROMode= " << mROMode << " to GRPUpdater"; pc.outputs().snapshot(Output{ "HMP", "ROMode", 0, Lifetime::Timeframe }, mROMode); @@ -150,6 +152,7 @@ o2::framework::DataProcessorSpec getHMPIDDigitizerSpec(int channel) Inputs{ InputSpec{ "collisioncontext", "SIM", "COLLISIONCONTEXT", static_cast(channel), Lifetime::Timeframe } }, Outputs{ OutputSpec{ "HMP", "DIGITS", 0, Lifetime::Timeframe }, + OutputSpec{ "HMP", "DIGITLBL", 0, Lifetime::Timeframe }, OutputSpec{ "HMP", "ROMode", 0, Lifetime::Timeframe } }, AlgorithmSpec{ adaptFromTask() },