AliceO2Group · sawenzel · Jan 23, 2019 · Jan 21, 2019 · Jan 22, 2019 · Jan 23, 2019
@@ -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,27 +93,27 @@ 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());
 
     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

@@ -9,8 +9,43 @@
 // or submit itself to any jurisdiction.
 
 #include "HMPIDBase/Digit.h"
-#include <iostream>
+#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;
+}
@@ -13,6 +13,8 @@
 
 #include "HMPIDBase/Digit.h"
 #include "HMPIDSimulation/Detector.h" // for the hit
+#include "SimulationDataFormat/MCTruthContainer.h"
+#include "SimulationDataFormat/MCCompLabel.h"
 #include <vector>
 
 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<o2::MCCompLabel>* labels)
+  {
+    mRegisteredLabelContainer = labels;
+  }
+
   // this will process hits and fill the digit vector with digits which are finalized
   void process(std::vector<o2::hmpid::HitType> const&, std::vector<o2::hmpid::Digit>& digit);
 
  private:
+  void zeroSuppress(std::vector<o2::hmpid::Digit> const& digits, std::vector<o2::hmpid::Digit>& newdigits,
+                    o2::dataformats::MCTruthContainer<o2::MCCompLabel> const& labels,
+                    o2::dataformats::MCTruthContainer<o2::MCCompLabel>* 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<o2::hmpid::Digit> mSummable;
   std::vector<o2::hmpid::Digit> mFinal;
 
-  // other stuff needed for digitizaton
+  std::vector<o2::hmpid::Digit> mDigits; // internal store for digits
+
+  //static constexpr int HMPID_NUMBEROFPADS = 161280;
+  //std::array<short, HMPID_NUMBEROFPADS> mIndexForPad = { -1 }; //! mapping of pad to digit index
+
+  std::map<int, short> mIndexForPad; //! logarithmic mapping of pad to digit index
+
+  std::vector<int> mInvolvedPads; //! list of pads where digits created
+
+  // other stuff needed for digitization
+  o2::dataformats::MCTruthContainer<o2::MCCompLabel> mTmpLabelContainer;                   // temp label container as workspace
+  o2::dataformats::MCTruthContainer<o2::MCCompLabel>* mRegisteredLabelContainer = nullptr; // label container to be filled
 
   ClassDefNV(HMPIDDigitizer, 1);
 };

@@ -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<o2::hmpid::Digit> const& digits, std::vector<o2::hmpid::Digit>& newdigits,
+                                  o2::dataformats::MCTruthContainer<o2::MCCompLabel> const& labels,
+                                  o2::dataformats::MCTruthContainer<o2::MCCompLabel>* 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<o2::hmpid::HitType> const& hits, std::vector<o2::hmpid::Digit>& 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<int, 9> 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);
 }
@@ -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