AliceO2Group · sawenzel · Aug 2, 2019 · Jul 31, 2019 · Jul 31, 2019 · Jul 31, 2019
@@ -83,13 +83,13 @@ class Digit : public DigitBase
     if (idx < kMaxLabels) {
       return mLabels[idx];
     } else {
-      return Label();
+      return -1;
     }
   }
   /// \brief Proportion of charge deposited by particle idx
   /// \param idx index in a list of a particles, max length kMaxLabels
   /// \return Proportion of energy from this particle.
-  double getLabelEProp(int idx) const
+  Label getLabelEProp(int idx) const
   {
     if (idx < kMaxLabels) {
       return mEProp[idx];

@@ -30,7 +30,6 @@ class Digit : public DigitBase
   using Label = o2::MCCompLabel;
 
  public:
-  static constexpr int kMaxLabels = 3; // Maximal number of MC labels associated with digit
   static constexpr int kTimeGate = 25; // Time in ns between digits to be added as one signal.
                                        // Should it be readout time (6000 ns???): to be tested
 
@@ -44,7 +43,7 @@ class Digit : public DigitBase
   /// \brief Digit constructor from Hit
   /// \param PHOS Hit
   /// \return constructed Digit
-  Digit(Hit hit);
+  Digit(Hit hit, int label);
 
   ~Digit() = default; // override
 
@@ -75,57 +74,36 @@ class Digit : public DigitBase
   void setAbsId(Int_t cellId) { mAbsId = cellId; }
 
   /// \brief Energy deposited in a cell
-  Double_t getAmplitude() const { return mAmplitude; }
-  void setAmplitude(Double_t amplitude) { mAmplitude = amplitude; }
+  float getAmplitude() const { return mAmplitude; }
+  void setAmplitude(float amplitude) { mAmplitude = amplitude; }
 
   /// \brief time measured in digit w.r.t. photon to PHOS arrival
-  Double_t getTime() const { return mTime; }
-  void setTime(Double_t time) { mTime = time; }
+  float getTime() const { return mTime; }
+  void setTime(float time) { mTime = time; }
 
   /// \brief Checks if this digit is produced in High Gain or Low Gain channels
-  Double_t isHighGain() const { return mIsHighGain; }
+  bool isHighGain() const { return mIsHighGain; }
   void setHighGain(Bool_t isHG) { mIsHighGain = isHG; }
 
-  /// \brief Label of a particle made energy deposition
-  /// \param idx index in a list of a particles, max length kMaxLabels
-  /// \return lable of a particle. Lables are sorted according to energy deposited by each of them
-  Label getLabel(Int_t idx) const
-  {
-    if (idx < kMaxLabels)
-      return mLabels[idx];
-    else
-      return Label();
-  }
-  /// \brief Proportion of energy deposited by particle idx
-  /// \param idx index in a list of a particles, max length kMaxLabels
-  /// \return Proportion of energy from this particle.
-  double getLabelEProp(Int_t idx) const
-  {
-    if (idx < kMaxLabels)
-      return mEProp[idx];
-    else
-      return 0.;
-  }
-  /// \brief Number of particles assosiated with this digit
-  int getNLabels() const { return mNlabels; }
+  /// \brief index of entry in MCLabels array
+  /// \return ndex of entry in MCLabels array
+  int getLabel() const { return mLabel; }
 
   void PrintStream(std::ostream& stream) const;
 
  private:
   // friend class boost::serialization::access;
 
-  int mAbsId;                ///< cell index (absolute cell ID)
-  double mAmplitude;         ///< Amplitude
-  double mTime;              ///< Time
-  int mNlabels;              ///< Number of actual labels in this digit
-  Label mLabels[kMaxLabels]; ///< Particle labels associated to this digit
-  double mEProp[kMaxLabels]; ///< Proportion of total energy deposited by given primary
-  bool mIsHighGain;          ///< High Gain or Low Gain channel (for calibration)
+  int mAbsId;       ///< cell index (absolute cell ID)
+  int mLabel;       ///< Index of the corresponding entry/entries in the MC label array
+  float mAmplitude; ///< Amplitude
+  float mTime;      ///< Time
+  bool mIsHighGain; ///< High Gain or Low Gain channel (for calibration)
 
   ClassDefNV(Digit, 1);
 };
 
 std::ostream& operator<<(std::ostream& stream, const Digit& dig);
-} // namespace PHOS
+} // namespace phos
 } // namespace o2
 #endif
@@ -15,42 +15,18 @@ using namespace o2::phos;
 
 ClassImp(Digit);
 
-constexpr int Digit::kMaxLabels;
-
 Digit::Digit(Int_t absId, Double_t amplitude, Double_t time, Int_t label)
-  : DigitBase(time), mAbsId(absId), mAmplitude(amplitude), mTime(time), mNlabels(0)
+  : DigitBase(time), mAbsId(absId), mAmplitude(amplitude), mTime(time), mLabel(label)
 {
-  if (label >= 0) {
-    mLabels[0] = label; // sofar there is no lables, no need to to sort
-    mEProp[0] = 1.;
-    mNlabels = 1;
-  }
 }
-Digit::Digit(Hit hit) : mAbsId(hit.GetDetectorID()), mAmplitude(hit.GetEnergyLoss()), mTime(hit.GetTime()), mNlabels(0)
+Digit::Digit(Hit hit, int label) : mAbsId(hit.GetDetectorID()), mAmplitude(hit.GetEnergyLoss()), mTime(hit.GetTime()), mLabel(label)
 {
-  mLabels[0] = hit.GetTrackID(); // so far there is no lables, no need to to sort
-  mEProp[0] = 1.;
-  for (Int_t i = 1; i < kMaxLabels; i++) {
-    mLabels[i] = -1;
-    mEProp[i] = 0.;
-  }
 }
 void Digit::FillFromHit(Hit hit)
 {
   mAbsId = hit.GetDetectorID();
   mAmplitude = hit.GetEnergyLoss();
   mTime = hit.GetTime();
-  if (hit.GetTrackID() >= 0) {
-    mLabels[0] = hit.GetTrackID(); // so far there is no lables, no need to to sort
-    mEProp[0] = 1.;
-    mNlabels = 1;
-  } else {
-    mNlabels = 0;
-  }
-  for (Int_t i = 1; i < kMaxLabels; i++) {
-    mLabels[i] = 0;
-    mEProp[i] = 0.;
-  }
 }
 
 bool Digit::operator<(const Digit& other) const
@@ -79,64 +55,7 @@ bool Digit::canAdd(const Digit other) const
 Digit& Digit::operator+=(const Digit& other)
 {
 
-  // Adds the amplitude of digits and completes the list of primary particles
-  double scaleThis = mAmplitude / (mAmplitude + other.mAmplitude);
-  double scaleOther = other.mAmplitude / (mAmplitude + other.mAmplitude);
-  for (Int_t i = 0; i < mNlabels; i++) {
-    mEProp[i] *= scaleThis;
-  }
-  if (other.mNlabels > 0) {
-    // copy and scale EProp of other digit
-    double otherEProp[kMaxLabels];
-    for (Int_t i = 0; i < other.mNlabels; i++) {
-      otherEProp[i] = scaleOther * other.mEProp[i];
-    }
-    double tmpEProp[kMaxLabels];
-    Label tmpLabels[kMaxLabels];
-
-    // Now find largest Energy Proportion
-    int i1 = 0, i2 = 0, i = 0;
-    while (i < kMaxLabels) {
-      if (i1 >= mNlabels) {
-        while (i2 < other.mNlabels) {
-          tmpEProp[i] = otherEProp[i2];
-          tmpLabels[i] = other.mLabels[i2];
-          i++;
-          i2++;
-        }
-        break;
-      }
-      if (i2 >= other.mNlabels) {
-        while (i1 < mNlabels) {
-          tmpEProp[i] = mEProp[i1];
-          tmpLabels[i] = mLabels[i1];
-          i++;
-          i1++;
-        }
-        break;
-      }
-
-      if (mEProp[i1] > otherEProp[i2]) {
-        tmpEProp[i] = mEProp[i1];
-        tmpLabels[i] = mLabels[i1];
-        i++;
-        i1++;
-      } else {
-        tmpEProp[i] = otherEProp[i2];
-        tmpLabels[i] = other.mLabels[i2];
-        i++;
-        i2++;
-      }
-    }
-    // Copy to current digit
-    for (Int_t ii = 0; ii < i; ii++) {
-      mEProp[ii] = tmpEProp[ii];
-      mLabels[ii] = tmpLabels[ii];
-    }
-  }
-
-  mNlabels = std::min(kMaxLabels, other.mNlabels + mNlabels);
-
+  // Adds the amplitude of digits
   // TODO: What about time? Should we assign time of more energetic digit? More complicated treatment?
   if (mAmplitude < other.mAmplitude) {
     mTime = other.mTime;

@@ -326,26 +326,13 @@ void Cluster::EvalPrimaries(const std::vector<Digit>* digits)
       continue;
     }
 
-    for (int iPrim = 0; iPrim < (*foundIt).getNLabels(); iPrim++) {
-      Label lab = (*foundIt).getLabel(iPrim);
-      double eProp = (*foundIt).getLabelEProp(iPrim);
-      // look through list of primaries, add proportion of energy if already there, else create new entry
-      int primListSize = mLabels.size();
-      bool found = false;
-      for (int inPrimList = 0; inPrimList < primListSize; inPrimList++) {
-        if (lab == mLabels[inPrimList]) { // already exist, increase proportin of energy
-          mLabelsEProp[inPrimList] += eProp * foundIt->getAmplitude() / mFullEnergy;
-          found = true;
-          break;
-        }
-      }
-      if (!found) {
-        mLabelsEProp.push_back(lab);
-        mLabelsEProp.push_back(eProp * ((*foundIt).getAmplitude()) / mFullEnergy);
-      }
-    }
+    //int lab = (*foundIt).getLabel();  //index of entry in MCLabels array
+    //Add Labels to list of primaries
+    //....
+    //TODO!!!!
   }
   // Vectors will not be modified any more.
+  //TODO: sort and add labels
   mLabels.shrink_to_fit();
   mLabelsEProp.shrink_to_fit();
 }

@@ -10,11 +10,13 @@
 
 o2_add_library(PHOSSimulation
                SOURCES src/Detector.cxx src/GeometryParams.cxx src/Digitizer.cxx
-                       src/DigitizerTask.cxx
+                       src/DigitizerTask.cxx src/PHOSSimParams.cxx
                PUBLIC_LINK_LIBRARIES O2::DetectorsBase O2::PHOSBase)
 
 o2_target_root_dictionary(PHOSSimulation
                           HEADERS include/PHOSSimulation/Detector.h
                                   include/PHOSSimulation/GeometryParams.h
                                   include/PHOSSimulation/Digitizer.h
-                                  include/PHOSSimulation/DigitizerTask.h)
+                                  include/PHOSSimulation/DigitizerTask.h
+                                  include/PHOSSimulation/PHOSSimParams.h
+                                  include/PHOSSimulation/MCLabel.h)
@@ -193,8 +193,8 @@ class Detector : public o2::base::DetImpl<Detector>
   friend class o2::base::DetImpl;
   ClassDefOverride(Detector, 1)
 };
-}
-}
+} // namespace phos
+} // namespace o2
 
 #ifdef USESHM
 namespace o2

@@ -14,6 +14,8 @@
 #include "PHOSBase/Digit.h"
 #include "PHOSBase/Geometry.h"
 #include "PHOSBase/Hit.h"
+#include "PHOSSimulation/MCLabel.h"
+#include "SimulationDataFormat/MCTruthContainer.h"
 
 namespace o2
 {
@@ -31,7 +33,7 @@ class Digitizer : public TObject
   void finish();
 
   /// Steer conversion of hits to digits
-  void process(const std::vector<Hit>& hits, std::vector<Digit>& digits);
+  void process(const std::vector<Hit>& hits, std::vector<Digit>& digits, o2::dataformats::MCTruthContainer<o2::phos::MCLabel>& labels);
 
   void setEventTime(double t);
   double getEventTime() const { return mEventTime; }
@@ -51,41 +53,24 @@ class Digitizer : public TObject
   void setCoeffToNanoSecond(double c) { mCoeffToNanoSecond = c; }
 
  protected:
-  Double_t NonLinearity(Double_t e);
-  Double_t DigitizeEnergy(Double_t e);
-  Double_t Decalibrate(Double_t e);
-  Double_t TimeResolution(Double_t time, Double_t e);
-  Double_t SimulateNoiseEnergy();
-  Double_t SimulateNoiseTime();
+  double NonLinearity(double e);
+  double DigitizeEnergy(double e);
+  double Decalibrate(double e);
+  double TimeResolution(double time, double e);
+  double SimulateNoiseEnergy();
+  double SimulateNoiseTime();
 
  private:
   const Geometry* mGeometry = nullptr; //!  PHOS geometry
   double mEventTime = 0;               ///< global event time
   double mCoeffToNanoSecond = 1.0;     ///< coefficient to convert event time (Fair) to ns
   bool mContinuous = false;            ///< flag for continuous simulation
-  UInt_t mROFrameMin = 0;              ///< lowest RO frame of current digits
-  UInt_t mROFrameMax = 0;              ///< highest RO frame of current digits
+  uint mROFrameMin = 0;                ///< lowest RO frame of current digits
+  uint mROFrameMax = 0;                ///< highest RO frame of current digits
   int mCurrSrcID = 0;                  ///< current MC source from the manager
   int mCurrEvID = 0;                   ///< current event ID from the manager
-  bool mApplyNonLinearity = true;      ///< if Non-linearity will be applied
-  bool mApplyDigitization = true;      ///< if energy digitization should be applied
-  bool mApplyDecalibration = false;    ///< if de-calibration should be applied
-  bool mApplyTimeResolution = true;    ///< if Hit time should be smeared
-  double mZSthreshold = 0.005;         ///< Zero Suppression threshold
-  double maNL = 0.04;                  ///< Parameter a for Non-Linearity
-  double mbNL = 0.2;                   ///< Parameter b for Non-Linearity
-  double mcNL = 1.;                    ///< Parameter c for Non-Linearity
-  double mADCWidth = 0.005;            ///< Widht of ADC channel used for energy digitization
-  double mTimeResolutionA = 2.;        ///< Time resolution parameter A (in ns)
-  double mTimeResolutionB = 2.;        ///< Time resolution parameter B (in ns/GeV)
-  double mTimeResThreshold = 0.5;      ///< threshold for time resolution calculation (in GeV)
-  double mAPDNoise = 0.005;            ///< Electronics (and APD) noise (in GeV)
-  double mMinNoiseTime = -200.;        ///< minimum time in noise channels (in ns)
-  double mMaxNoiseTime = 2000.;        ///< minimum time in noise channels (in ns)
 
-  //  std::unordered_map<Int_t, std::deque<Digit>> mDigits; ///< used to sort digits by tower
-
-  ClassDefOverride(Digitizer, 1);
+  ClassDefOverride(Digitizer, 2);
 };
 } // namespace phos
 } // namespace o2

@@ -38,14 +38,15 @@ class DigitizerTask : public FairTask
   double getFairTimeUnitInNS() const { return mFairTimeUnitInNS; }
 
  private:
-  double mFairTimeUnitInNS = 1;                 ///< Fair time unit in ns
-  Int_t mSourceID = 0;                          ///< current source
-  Int_t mEventID = 0;                           ///< current event id from the source
-  Digitizer mDigitizer;                         ///< Digitizer
-  const std::vector<Hit>* mHitsArray = nullptr; ///< Array of MC hits
-  std::vector<Digit>* mDigitsArray = nullptr;   ///< Array of digits
-
-  ClassDefOverride(DigitizerTask, 1);
+  double mFairTimeUnitInNS = 1;                                            ///< Fair time unit in ns
+  Int_t mSourceID = 0;                                                     ///< current source
+  Int_t mEventID = 0;                                                      ///< current event id from the source
+  Digitizer mDigitizer;                                                    ///< Digitizer
+  const std::vector<Hit>* mHitsArray = nullptr;                            ///< Array of MC hits
+  std::vector<Digit>* mDigitsArray = nullptr;                              ///< Array of digits
+  o2::dataformats::MCTruthContainer<o2::phos::MCLabel>* mLabels = nullptr; ///< Array of digit labels
+
+  ClassDefOverride(DigitizerTask, 2);
 };
 } // namespace phos
 } // namespace o2