T0 digitzation as a part of FIT

AliceO2Group · sawenzel · Apr 1, 2019 · Jan 18, 2019 · Jan 22, 2019 · Jan 23, 2019
commit 9abcec91014a43e54464f31483334ab6a3fd2fd8
@@ -26,7 +26,6 @@ class FairVolume;
 class TGeoVolume;
 class TGraph;
 
-
 namespace o2
 {
 namespace t0

@@ -19,6 +19,7 @@ inline o2::fit::DigitizationParameters T0DigitizationParameters()
   result.mTime_trg_gate = 4.;                                 // ns
   result.mAmpThreshold = 100;                                 // number of photoelectrons
   result.mTimeDiffAC = (Geometry::ZdetA - Geometry::ZdetC) * TMath::C();
+  result.mIsT0 = true;
 
   return result;
 }

@@ -59,12 +59,12 @@ class DigitizerTask : public FairTask
   Digitizer mDigitizer;                                     ///< Digitizer
   const std::vector<o2::t0::HitType>* mHitsArray = nullptr; ///< Array of MC hits
 
-  TBranch* mQEDBranch = nullptr;                               //! optional special branch of hits from QED collitions
+  TBranch* mQEDBranch = nullptr;                                //! optional special branch of hits from QED collitions
   const std::vector<o2::fit::HitType>* mHitsArrayQED = nullptr; //! array of MC hits from ED
-  float mQEDEntryTimeBinNS = 0.f;                              ///< every entry in the QED branch integrates QED for so many nanosec.
-  double mLastQEDTimeNS = 0.;                                  ///< center of the time-bin of last added QED bg slot (entry of mQEDBranch)
-  int mLastQEDEntry = -1;                                      ///< last used QED entry
-  UChar_t mQEDSourceID = 0;                                    ///< MC ID source of the QED (stored in the labels)
+  float mQEDEntryTimeBinNS = 0.f;                               ///< every entry in the QED branch integrates QED for so many nanosec.
+  double mLastQEDTimeNS = 0.;                                   ///< center of the time-bin of last added QED bg slot (entry of mQEDBranch)
+  int mLastQEDEntry = -1;                                       ///< last used QED entry
+  UChar_t mQEDSourceID = 0;                                     ///< MC ID source of the QED (stored in the labels)
 
   o2::fit::Digit* mEventDigit = nullptr;
   o2::dataformats::MCTruthContainer<o2::fit::MCLabel> mMCTruthArray;                      //! Labels containter

@@ -34,7 +34,7 @@ using o2::t0::Geometry;
 ClassImp(Detector);
 
 Detector::Detector(Bool_t Active)
-  : o2::Base::DetImpl<Detector> ("T0", Active), mIdSens1(0), mPMTeff(nullptr), mHits(o2::utils::createSimVector<o2::fit::HitType>())
+  : o2::Base::DetImpl<Detector>("T0", Active), mIdSens1(0), mPMTeff(nullptr), mHits(o2::utils::createSimVector<o2::fit::HitType>())
 
 {
   // Gegeo  = GetGeometry() ;
@@ -43,7 +43,7 @@ Detector::Detector(Bool_t Active)
 }
 
 Detector::Detector(const Detector& rhs)
-  : o2::Base::DetImpl<Detector> (rhs), mIdSens1(rhs.mIdSens1), mPMTeff(rhs.mPMTeff), mHits(o2::utils::createSimVector<o2::fit::HitType>())
+  : o2::Base::DetImpl<Detector>(rhs), mIdSens1(rhs.mIdSens1), mPMTeff(rhs.mPMTeff), mHits(o2::utils::createSimVector<o2::fit::HitType>())
 {
 }
 

@@ -24,7 +24,8 @@ struct ChannelData {
   Int_t ChId;       //channel Id
   Double_t CFDTime; //time in ns, 0 at lhc clk center
   Double_t QTCAmpl; // Amplitude in mips
-  ClassDefNV(ChannelData, 1);
+  int numberOfParticles;
+  ClassDefNV(ChannelData, 2);
 };
 
 /// \class Digit
@@ -71,6 +72,7 @@ class Digit : public DigitBase
   }
 
   const std::vector<ChannelData>& getChDgData() const { return mChDgDataArr; }
+  std::vector<ChannelData>& getChDgData() { return mChDgDataArr; }
   void setChDgData(const std::vector<ChannelData>& ChDgDataArr) { mChDgDataArr = ChDgDataArr; }
   void setChDgData(std::vector<ChannelData>&& ChDgDataArr) { mChDgDataArr = std::move(ChDgDataArr); }
 

@@ -16,6 +16,7 @@ struct DigitizationParameters {
   Float_t mTime_trg_gate; // ns
   Int_t mAmpThreshold;    // number of photoelectrons
   Float_t mTimeDiffAC;
+  bool mIsT0; //amplitude T0(true) or V0 (false)
 };
 } // namespace o2::fit
 #endif
@@ -30,6 +30,7 @@ class Digitizer
 
   //void process(const std::vector<HitType>* hits, std::vector<Digit>* digits);
   void process(const std::vector<o2::fit::HitType>* hits, Digit* digit);
+  void computeAverage(Digit& digit);
 
   void initParameters();
   // void printParameters();
@@ -58,14 +59,6 @@ class Digitizer
   Double_t mEventTime; // timestamp
 
   DigitizationParameters parameters;
-  /*
-  Float_t mBC_clk_center = 12.5;                             // clk center
-  Int_t mMCPs = (Geometry::NCellsA + Geometry::NCellsC) * 4; //number of MCPs
-  Float_t mCFD_trsh_mip = 0.4;                               // = 4[mV] / 10[mV/mip]
-  Float_t mTime_trg_gate = 4.;                               // ns
-  Int_t mAmpThreshold = 100;                                 // number of photoelectrons
-  Float_t mTimeDiffAC = (Geometry::ZdetA - Geometry::ZdetC) * TMath::C();
-  */
 
   o2::dataformats::MCTruthContainer<o2::fit::MCLabel>* mMCLabels = nullptr;
 

@@ -14,16 +14,16 @@ class HitType : public o2::BasicXYZEHit<float>
  public:
   using BasicXYZEHit<float>::BasicXYZEHit;
 };
-} // namespace t0
+} // namespace fit
 } // namespace o2
 
 #endif
 
 #ifdef USESHM
 namespace std
 {
-  template <>
-  class allocator<o2::fit::HitType> : public o2::utils::ShmAllocator<o2::fit::HitType>
+template <>
+class allocator<o2::fit::HitType> : public o2::utils::ShmAllocator<o2::fit::HitType>
 {
 };
 } // namespace std

@@ -28,10 +28,9 @@ void Digitizer::process(const std::vector<o2::fit::HitType>* hits, Digit* digit)
 {
   //parameters constants TO DO: move to class
 
-  constexpr Float_t C_side_cable_cmps = 2.8; //ns
-  constexpr Float_t A_side_cable_cmps = 11.; //ns
-  constexpr Float_t signal_width = 5.;       // time gate for signal, ns
-  constexpr Float_t nPe_in_mip = 250.;       // n ph. e. in one mip
+  constexpr Float_t C_side_cable_cmps = 2.877; //ns
+  constexpr Float_t A_side_cable_cmps = 11.08; //ns
+  constexpr Float_t signal_width = 5.;         // time gate for signal, ns
 
   Int_t nlbl = 0; //number of MCtrues
 
@@ -40,11 +39,13 @@ void Digitizer::process(const std::vector<o2::fit::HitType>* hits, Digit* digit)
   digit->setOrbit(mOrbit);
 
   //Calculating signal time, amplitude in mean_time +- time_gate --------------
-  Double_t cfd[300] = {};
-  Float_t amp[300] = {};
-  Int_t ch_signal_nPe[300] = {};
-  Double_t ch_signal_MIP[300] = {};
-  Double_t ch_signal_time[300] = {};
+  std::vector<ChannelData>& channel_data = digit->getChDgData();
+  if (channel_data.size() == 0) {
+    channel_data.reserve(parameters.mMCPs);
+    for (int i = 0; i < parameters.mMCPs; ++i)
+      channel_data.emplace_back(ChannelData{ i, 0, 0, 0 });
+  }
+  assert(digit->getChDgData().size() == parameters.mMCPs);
   for (auto& hit : *hits) {
     Int_t hit_ch = hit.GetDetectorID();
     Double_t hit_time = hit.GetTime();
@@ -56,9 +57,10 @@ void Digitizer::process(const std::vector<o2::fit::HitType>* hits, Digit* digit)
 
     Double_t hit_time_corr = hit_time - time_compensate /* + mBC_clk_center + mEventTime*/;
 
-    if (/*is_time_in_gate &&*/ is_hit_in_signal_gate) {
-      ch_signal_nPe[hit_ch]++;
-      ch_signal_time[hit_ch] += hit_time_corr;
+    if (is_hit_in_signal_gate) {
+      channel_data[hit_ch].numberOfParticles++;
+      channel_data[hit_ch].QTCAmpl += hit.GetEnergyLoss();
+      channel_data[hit_ch].CFDTime += hit_time_corr;
     }
 
     //charge particles in MCLabel
@@ -72,37 +74,24 @@ void Digitizer::process(const std::vector<o2::fit::HitType>* hits, Digit* digit)
       }
     }
   }
+}
 
-  // sum  different sources
-  std::vector<ChannelData> mChDgDataArr;
-  for (const auto& d : digit->getChDgData()) {
-    Int_t mcp = d.ChId;
-    cfd[mcp] = d.CFDTime;
-    amp[mcp] = d.QTCAmpl;
-  }
-
-  for (Int_t ch_iter = 0; ch_iter < parameters.mMCPs; ch_iter++) {
-    if (ch_signal_nPe[ch_iter] != 0) {
-      ch_signal_MIP[ch_iter] = amp[ch_iter] + ch_signal_nPe[ch_iter] / nPe_in_mip;
-      if (cfd[ch_iter] > 0) {
-        cfd[ch_iter] = cfd[ch_iter] - parameters.mBC_clk_center - mEventTime;
-        ch_signal_time[ch_iter] = ((cfd[ch_iter] + ch_signal_time[ch_iter] / (float)ch_signal_nPe[ch_iter]) / 2.) + parameters.mBC_clk_center + mEventTime;
-      } else
-        ch_signal_time[ch_iter] = (ch_signal_time[ch_iter] / (float)ch_signal_nPe[ch_iter]) + parameters.mBC_clk_center + mEventTime;
-
-      if (ch_signal_MIP[ch_iter] > parameters.mCFD_trsh_mip) {
-        LOG(DEBUG) << ch_iter << " : "
-                   << " : " << ch_signal_time[ch_iter] - parameters.mBC_clk_center - mEventTime << " : "
-                   << ch_signal_MIP[ch_iter] << " : " << mEventTime << " cfd " << cfd[ch_iter] << FairLogger::endl;
-        mChDgDataArr.emplace_back(ChannelData{ ch_iter, ch_signal_time[ch_iter], ch_signal_MIP[ch_iter] });
-      } else {
-        ch_signal_MIP[ch_iter] = 0;
-        ch_signal_time[ch_iter] = 0;
-      }
-    }
+void Digitizer::computeAverage(Digit& digit)
+{
+  constexpr Float_t nPe_in_mip = 250.; // n ph. e. in one mip
+  auto& channel_data = digit.getChDgData();
+  for (auto& ch_data : channel_data) {
+    if (ch_data.numberOfParticles == 0)
+      continue;
+    ch_data.CFDTime /= ch_data.numberOfParticles;
+    if (parameters.mIsT0)
+      ch_data.QTCAmpl = ch_data.numberOfParticles / nPe_in_mip;
   }
-
-  digit->setChDgData(std::move(mChDgDataArr));
+  channel_data.erase(std::remove_if(channel_data.begin(), channel_data.end(),
+                                    [this](ChannelData const& ch_data) {
+                                      return ch_data.QTCAmpl < parameters.mCFD_trsh_mip;
+                                    }),
+                     channel_data.end());
 }
 
 //------------------------------------------------------------------------
@@ -112,11 +101,12 @@ void Digitizer::smearCFDtime(Digit* digit)
   std::vector<ChannelData> mChDgDataArr;
   for (const auto& d : digit->getChDgData()) {
     Int_t mcp = d.ChId;
-    Double_t cfd = d.CFDTime - parameters.mBC_clk_center - mEventTime;
+    Double_t cfd = d.CFDTime;
     Float_t amp = d.QTCAmpl;
+    int numpart = d.numberOfParticles;
     if (amp > parameters.mCFD_trsh_mip) {
       Double_t smeared_time = gRandom->Gaus(cfd, 0.050) + parameters.mBC_clk_center + mEventTime;
-      mChDgDataArr.emplace_back(ChannelData{ mcp, smeared_time, amp });
+      mChDgDataArr.emplace_back(ChannelData{ mcp, smeared_time, amp, numpart });
     }
   }
   digit->setChDgData(std::move(mChDgDataArr));
@@ -143,7 +133,7 @@ void Digitizer::setTriggers(Digit* digit)
   Float_t amp[300] = {};
   for (const auto& d : digit->getChDgData()) {
     Int_t mcp = d.ChId;
-    cfd[mcp] = d.CFDTime - parameters.mBC_clk_center - mEventTime;
+    cfd[mcp] = d.CFDTime;
     amp[mcp] = d.QTCAmpl;
     if (amp[mcp] < parameters.mCFD_trsh_mip)
       continue;
@@ -204,8 +194,9 @@ void Digitizer::initParameters()
   // murmur
 }
 //_______________________________________________________________________
-void Digitizer::init() {
-  std::cout<<" @@@ Digitizer::init "<<std::endl;
+void Digitizer::init()
+{
+  std::cout << " @@@ Digitizer::init " << std::endl;
 }
 
 //_______________________________________________________________________

@@ -25,7 +25,7 @@
 #include <memory> // for make_shared, make_unique, unique_ptr
 #include <vector>
 #include <string>
-#include <algorithm> 
+#include <algorithm>
 
 using namespace o2::framework;
 using SubSpecificationType = o2::framework::DataAllocator::SubSpecificationType;
@@ -44,11 +44,11 @@ class FITDPLDigitWriter
   void init(framework::InitContext& ic)
   {
     std::string detStrL = mID.getName();
-     std::transform(detStrL.begin(), detStrL.end(), detStrL.begin(), ::tolower);
+    std::transform(detStrL.begin(), detStrL.end(), detStrL.begin(), ::tolower);
 
     auto filename = ic.options().get<std::string>((detStrL + "-digit-outfile").c_str());
     auto treename = ic.options().get<std::string>("treename");
-    std::cout<<" @@@ FITDPLDigitWriter " <<detStrL<<" file "<<filename<<" tree "<<treename <<std::endl;
+    std::cout << " @@@ FITDPLDigitWriter " << detStrL << " file " << filename << " tree " << treename << std::endl;
 
     mOutFile = std::make_unique<TFile>(filename.c_str(), "RECREATE");
     if (!mOutFile || mOutFile->IsZombie()) {
@@ -154,6 +154,5 @@ DataProcessorSpec getT0DigitWriterSpec()
   };
 }
 
-
 } // end namespace fit
 } // end namespace o2
@@ -21,7 +21,7 @@ namespace fit
 o2::framework::DataProcessorSpec getT0DigitWriterSpec();
 //o2::framework::DataProcessorSpec getMFTDigitWriterSpec();
 
-} // end namespace itsmft
+} // namespace fit
 } // end namespace o2
- 
+
 #endif /* STEER_DIGITIZERWORKFLOW_ITSMFTDIGITWRITER_H_ */