From 10cbc2698a45b9b321d5f96f38d1f18fd7adea21 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Mon, 21 Jan 2019 16:52:22 +0100 Subject: [PATCH 1/4] Improvded HMPID Digitizer --- .../HMPID/base/include/HMPIDBase/Digit.h | 59 +++++++++++-------- Detectors/HMPID/base/src/Digit.cxx | 37 +++++++++++- .../include/HMPIDSimulation/HMPIDDigitizer.h | 7 ++- .../HMPID/simulation/src/HMPIDDigitizer.cxx | 53 ++++++++++++++++- 4 files changed, 128 insertions(+), 28 deletions(-) diff --git a/Detectors/HMPID/base/include/HMPIDBase/Digit.h b/Detectors/HMPID/base/include/HMPIDBase/Digit.h index a19aeb33b2f7e..556b413e8073c 100644 --- a/Detectors/HMPID/base/include/HMPIDBase/Digit.h +++ b/Detectors/HMPID/base/include/HMPIDBase/Digit.h @@ -32,47 +32,56 @@ class Digit : public DigitBase float getCharge() const { return mQ; } int getPadID() const { return 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 +90,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,9 +105,9 @@ 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); 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..170d1dbe4561e 100644 --- a/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h +++ b/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h @@ -39,7 +39,12 @@ class HMPIDDigitizer std::vector mSummable; std::vector mFinal; - // other stuff needed for digitizaton + std::vector mDigits; + static constexpr int HMPID_NUMBEROFPADS = 161280; + std::array mIndexForPad = { -1 }; //! mapping of pad to digit index + std::vector mInvolvedPads; //! list of pads where digits created + + // other stuff needed for digitization ClassDefNV(HMPIDDigitizer, 1); }; diff --git a/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx b/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx index d4d785d724ec1..dc0e1a862883e 100644 --- a/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx +++ b/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx @@ -21,7 +21,58 @@ void HMPIDDigitizer::process(std::vector const& hits, std::v // 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. + + // for (auto& hit : hits) { + // digits.emplace_back(hit); + // } + + // clear lookup structures + for (auto& pad : mInvolvedPads) { + mIndexForPad[pad] = -1; + } + 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); + LOG(INFO) << "CHAMBER " << chamber; + LOG(INFO) << "PC " << pc; + LOG(INFO) << "PX " << px; + LOG(INFO) << "PY " << py; + + 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); + LOG(INFO) << "ADDING PAD " << allpads[counter]; + counter++; + } + } + + LOG(INFO) << "DIGIT ON MAINPAD " << Param::Abs(chamber, pc, px, py) << " TOTAL CHARGE " << totalQ; + + for (auto& pad : allpads) { + 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); + } else { + // create digit ... and register + mDigits.emplace_back(pad, totalQ * fraction); + mIndexForPad[pad] = mDigits.size() - 1; + mInvolvedPads.emplace_back(pad); + } + } } } From ba0a3882f370e3e689f892fe987a78766e0cf82c Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Tue, 22 Jan 2019 13:17:18 +0100 Subject: [PATCH 2/4] Digitizer: Threshold; Correct index mapping --- .../HMPID/base/include/HMPIDBase/Digit.h | 5 +- .../include/HMPIDSimulation/HMPIDDigitizer.h | 15 ++++-- .../HMPID/simulation/src/HMPIDDigitizer.cxx | 46 ++++++++++++++----- Steer/DigitizerWorkflow/CMakeLists.txt | 2 +- .../src/HMPIDDigitWriterSpec.h | 23 ++++++++-- 5 files changed, 72 insertions(+), 19 deletions(-) diff --git a/Detectors/HMPID/base/include/HMPIDBase/Digit.h b/Detectors/HMPID/base/include/HMPIDBase/Digit.h index 556b413e8073c..381f2556b273b 100644 --- a/Detectors/HMPID/base/include/HMPIDBase/Digit.h +++ b/Detectors/HMPID/base/include/HMPIDBase/Digit.h @@ -31,6 +31,9 @@ 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); } static void getPadAndTotalCharge(HitType const& hit, int& chamber, int& pc, int& px, int& py, float& totalcharge) { @@ -110,7 +113,7 @@ class Digit : public DigitBase return 4. * Digit::IntPartMathiX(localX, pad) * Digit::IntPartMathiY(localY, pad); } - ClassDefNV(Digit, 1); + ClassDefNV(Digit, 2); }; } // namespace hmpid diff --git a/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h b/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h index 170d1dbe4561e..e9a94c1ae7105 100644 --- a/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h +++ b/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h @@ -31,6 +31,11 @@ class HMPIDDigitizer void process(std::vector const&, std::vector& digit); private: + + void zeroSuppress(std::vector const& digits, std::vector& newdigits); + 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,9 +44,13 @@ class HMPIDDigitizer std::vector mSummable; std::vector mFinal; - std::vector mDigits; - static constexpr int HMPID_NUMBEROFPADS = 161280; - std::array mIndexForPad = { -1 }; //! mapping of pad to digit index + 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 diff --git a/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx b/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx index dc0e1a862883e..a7b2192459831 100644 --- a/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx +++ b/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx @@ -15,6 +15,20 @@ 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) { + for(auto& digit: digits) { + if(digit.getCharge() >= getThreshold(digit)) { + newdigits.push_back(digit); + } + } +} + // this will process hits and fill the digit vector with digits which are finalized void HMPIDDigitizer::process(std::vector const& hits, std::vector& digits) { @@ -27,9 +41,10 @@ void HMPIDDigitizer::process(std::vector const& hits, std::v // } // clear lookup structures - for (auto& pad : mInvolvedPads) { - mIndexForPad[pad] = -1; - } + //for (auto& pad : mInvolvedPads) { +// mIndexForPad[pad] = -1; +// } + mIndexForPad.clear(); mInvolvedPads.clear(); for (auto& hit : hits) { @@ -59,20 +74,29 @@ void HMPIDDigitizer::process(std::vector const& hits, std::v LOG(INFO) << "DIGIT ON MAINPAD " << Param::Abs(chamber, pc, px, py) << " TOTAL CHARGE " << totalQ; + for (auto& pad : allpads) { - auto index = mIndexForPad[pad]; + 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); + // digit exists ... reuse + auto& digit = mDigits[index]; + digit.addCharge(totalQ * fraction); } else { - // create digit ... and register - mDigits.emplace_back(pad, totalQ * fraction); - mIndexForPad[pad] = mDigits.size() - 1; - mInvolvedPads.emplace_back(pad); + // create digit ... and register + mDigits.emplace_back(pad, totalQ * fraction); + mIndexForPad[pad] = mDigits.size() - 1; + mInvolvedPads.emplace_back(pad); } } } + + // + zeroSuppress(mDigits, digits); } diff --git a/Steer/DigitizerWorkflow/CMakeLists.txt b/Steer/DigitizerWorkflow/CMakeLists.txt index 03d41f1c1b43b..4642e6f207ee9 100644 --- a/Steer/DigitizerWorkflow/CMakeLists.txt +++ b/Steer/DigitizerWorkflow/CMakeLists.txt @@ -37,7 +37,7 @@ O2_GENERATE_EXECUTABLE( src/EMCALDigitizerSpec.cxx src/EMCALDigitWriterSpec.cxx src/HMPIDDigitizerSpec.cxx - src/HMPIDDigitWriterSpec.cxx + # src/HMPIDDigitWriterSpec.cxx src/MCHDigitizerSpec.cxx src/TRDDigitizerSpec.cxx src/GRPUpdaterSpec.cxx diff --git a/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.h b/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.h index c03355d580160..1435d45f59b0e 100644 --- a/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.h +++ b/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.h @@ -8,17 +8,34 @@ // 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" 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{ "input", "HMP", "DIGITS" }, "HMPDigit" } + // add more branch definitions (for example Monte Carlo labels here) + )(); +} } // end namespace hmpid } // end namespace o2 From 192b0138dab008ce576167fccb7a3b6c754b2fae Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Wed, 23 Jan 2019 11:53:50 +0100 Subject: [PATCH 3/4] Implemented labels --- .../include/HMPIDSimulation/HMPIDDigitizer.h | 15 ++- .../HMPID/simulation/src/HMPIDDigitizer.cxx | 95 +++++++++++-------- .../src/HMPIDDigitWriterSpec.h | 7 +- .../src/HMPIDDigitizerSpec.cxx | 11 ++- 4 files changed, 78 insertions(+), 50 deletions(-) diff --git a/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h b/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h index e9a94c1ae7105..0e0666a3519ae 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,12 +29,21 @@ 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); - void zeroSuppress(std::vector const& digits, std::vector& newdigits); float getThreshold(o2::hmpid::Digit const&) const; // gives back threshold to apply for a certain digit // (using noise and other tables for pad) @@ -54,6 +65,8 @@ class HMPIDDigitizer 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 a7b2192459831..925efc4f1bbf9 100644 --- a/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx +++ b/Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx @@ -15,35 +15,34 @@ using namespace o2::hmpid; ClassImp(HMPIDDigitizer); -float HMPIDDigitizer::getThreshold(o2::hmpid::Digit const& digiti) const { +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) { - for(auto& digit: digits) { - if(digit.getCharge() >= getThreshold(digit)) { - newdigits.push_back(digit); - } - } +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. - - // for (auto& hit : hits) { - // digits.emplace_back(hit); - // } - - // clear lookup structures - //for (auto& pad : mInvolvedPads) { -// mIndexForPad[pad] = -1; -// } mIndexForPad.clear(); mInvolvedPads.clear(); @@ -52,10 +51,6 @@ void HMPIDDigitizer::process(std::vector const& hits, std::v float totalQ; // retrieves center pad and the total charge Digit::getPadAndTotalCharge(hit, chamber, pc, px, py, totalQ); - LOG(INFO) << "CHAMBER " << chamber; - LOG(INFO) << "PC " << pc; - LOG(INFO) << "PX " << px; - LOG(INFO) << "PY " << py; if (px < 0 || py < 0) { continue; @@ -67,36 +62,52 @@ void HMPIDDigitizer::process(std::vector const& hits, std::v for (int nx = -1; nx <= 1; ++nx) { for (int ny = -1; ny <= 1; ++ny) { allpads[counter] = Param::Abs(chamber, pc, px + nx, py + ny); - LOG(INFO) << "ADDING PAD " << allpads[counter]; counter++; } } - LOG(INFO) << "DIGIT ON MAINPAD " << Param::Abs(chamber, pc, px, py) << " TOTAL CHARGE " << totalQ; - - for (auto& pad : allpads) { - auto iter = mIndexForPad.find(pad); - int index = -1; - if (iter != mIndexForPad.end()) { - index = iter->second; - } + 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; + // LOG(INFO) << "FRACTION ON PAD " << pad << " IS " << fraction; if (index != -1) { - // digit exists ... reuse - auto& digit = mDigits[index]; - digit.addCharge(totalQ * fraction); + // 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); + // 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)); + } } } } - // - zeroSuppress(mDigits, digits); + // apply zero suppression (removes digits in the noise) + zeroSuppress(mDigits, digits, mTmpLabelContainer, mRegisteredLabelContainer); } diff --git a/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.h b/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.h index 1435d45f59b0e..4b26b7d60fe49 100644 --- a/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.h +++ b/Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.h @@ -15,6 +15,8 @@ #include "Utils/MakeRootTreeWriterSpec.h" #include "Framework/InputSpec.h" #include "HMPIDBase/Digit.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "SimulationDataFormat/MCCompLabel.h" namespace o2 { @@ -32,9 +34,8 @@ o2::framework::DataProcessorSpec getHMPIDDigitWriterSpec() "hmpiddigits.root", "o2sim", 1, - BranchDefinition>{ InputSpec{ "input", "HMP", "DIGITS" }, "HMPDigit" } - // add more branch definitions (for example Monte Carlo labels here) - )(); + BranchDefinition>{ InputSpec{ "digitinput", "HMP", "DIGITS" }, "HMPDigit" }, + BranchDefinition>{ InputSpec{ "labelinput", "HMP", "DIGITLBL" }, "HMPDigitLabels" })(); } } // end namespace hmpid 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() }, From ba1894beaf922316b53184319f5d92dfcda1da15 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Wed, 23 Jan 2019 11:55:57 +0100 Subject: [PATCH 4/4] cleanup --- .../include/HMPIDSimulation/HMPIDDigitizer.h | 4 +- Steer/DigitizerWorkflow/CMakeLists.txt | 1 - .../src/HMPIDDigitWriterSpec.cxx | 73 ------------------- 3 files changed, 2 insertions(+), 76 deletions(-) delete mode 100644 Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.cxx diff --git a/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h b/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h index 0e0666a3519ae..38b26e4870d40 100644 --- a/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h +++ b/Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h @@ -45,7 +45,7 @@ class HMPIDDigitizer 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) + // (using noise and other tables for pad) double mTime = 0.; int mEventID = 0; @@ -65,7 +65,7 @@ class HMPIDDigitizer 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 mTmpLabelContainer; // temp label container as workspace o2::dataformats::MCTruthContainer* mRegisteredLabelContainer = nullptr; // label container to be filled ClassDefNV(HMPIDDigitizer, 1); diff --git a/Steer/DigitizerWorkflow/CMakeLists.txt b/Steer/DigitizerWorkflow/CMakeLists.txt index 4642e6f207ee9..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