add sipm digitization

simulation/g4simulation/g4calo/RawTowerDigitizer.cc

@@ -1,35 +1,36 @@
 #include "RawTowerDigitizer.h"
 
-#include <calobase/RawTower.h>               // for RawTower
-#include <calobase/RawTowerv1.h>
+#include <calobase/RawTower.h>  // for RawTower
 #include <calobase/RawTowerContainer.h>
 #include <calobase/RawTowerDeadMap.h>
-#include <calobase/RawTowerDefs.h>           // for keytype
+#include <calobase/RawTowerDefs.h>  // for keytype
 #include <calobase/RawTowerGeom.h>
 #include <calobase/RawTowerGeomContainer.h>
+#include <calobase/RawTowerv1.h>
 
-#include <fun4all/Fun4AllBase.h>             // for Fun4AllBase::VERBOSITY_MORE
+#include <fun4all/Fun4AllBase.h>  // for Fun4AllBase::VERBOSITY_MORE
 #include <fun4all/Fun4AllReturnCodes.h>
-#include <fun4all/SubsysReco.h>              // for SubsysReco
+#include <fun4all/SubsysReco.h>  // for SubsysReco
 
 #include <phool/PHCompositeNode.h>
 #include <phool/PHIODataNode.h>
-#include <phool/PHNode.h>                    // for PHNode
+#include <phool/PHNode.h>  // for PHNode
 #include <phool/PHNodeIterator.h>
-#include <phool/PHObject.h>                  // for PHObject
+#include <phool/PHObject.h>  // for PHObject
 #include <phool/PHRandomSeed.h>
 #include <phool/getClass.h>
 
 #include <gsl/gsl_randist.h>
 #include <gsl/gsl_rng.h>
+#include <gsl/gsl_cdf.h>
 
 #include <cmath>
-#include <cstdlib>                          // for exit
-#include <exception>                         // for exception
+#include <cstdlib>    // for exit
+#include <exception>  // for exception
 #include <iostream>
 #include <map>
 #include <stdexcept>
-#include <utility>                           // for pair
+#include <utility>  // for pair
 
 using namespace std;
 
@@ -49,6 +50,7 @@ RawTowerDigitizer::RawTowerDigitizer(const std::string &name)
   , m_PedstalWidthADC(NAN)
   , m_ZeroSuppressionADC(0)  //default to apply no zero suppression
   , m_TowerType(-1)
+  , m_SiPMEffectivePixel(40000 * 4)  // sPHENIX EMCal default, 4x Hamamatsu S12572-015P MPPC [sPHENIX TDR]
 {
   m_RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
   m_Seed = PHRandomSeed();  // fixed seed handled in PHRandomSeed()
@@ -162,6 +164,11 @@ int RawTowerDigitizer::process_event(PHCompositeNode *topNode)
       // for photon digitization towers can be created if sim_tower is null pointer
       digi_tower = simple_photon_digitization(sim_tower);
     }
+    else if (m_DigiAlgorithm == kSiPM_photon_digitization)
+    {
+      // for photon digitization towers can be created if sim_tower is null pointer
+      digi_tower = simple_photon_digitization(sim_tower);
+    }
     else
     {
       cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
@@ -257,6 +264,72 @@ RawTowerDigitizer::simple_photon_digitization(RawTower *sim_tower)
   return digi_tower;
 }
 
+RawTower *
+RawTowerDigitizer::sipm_photon_digitization(RawTower *sim_tower)
+{
+  RawTower *digi_tower = nullptr;
+
+  double energy = 0;
+  if (sim_tower)
+  {
+    energy = sim_tower->get_energy();
+  }
+  const double photon_count_mean = energy * m_PhotonElecYieldVisibleGeV;
+  const double poission_param_per_pixel = photon_count_mean / m_SiPMEffectivePixel;
+  const double prob_activated_per_pixel = gsl_cdf_poisson_Q(0, poission_param_per_pixel);
+  const double active_pixel = gsl_ran_binomial(m_RandomGenerator , prob_activated_per_pixel, m_SiPMEffectivePixel);
+
+//  const int photon_count = gsl_ran_poisson(m_RandomGenerator, photon_count_mean);
+
+  const int signal_ADC = floor(active_pixel / m_PhotonElecADC);
+
+  const double pedstal = m_PedstalCentralADC + ((m_PedstalWidthADC > 0) ? gsl_ran_gaussian(m_RandomGenerator, m_PedstalWidthADC) : 0);
+  const int sum_ADC = signal_ADC + (int) pedstal;
+
+  if (sum_ADC > m_ZeroSuppressionADC)
+  {
+    // create new digitalizaed tower
+    if (sim_tower)
+    {
+      digi_tower = new RawTowerv1(*sim_tower);
+    }
+    else
+    {
+      digi_tower = new RawTowerv1();
+    }
+    digi_tower->set_energy((double) sum_ADC);
+  }
+
+  if (Verbosity() >= 2)
+  {
+    cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
+         << endl;
+
+    cout << "input: ";
+    if (sim_tower)
+    {
+      sim_tower->identify();
+    }
+    else
+    {
+      cout << "None" << endl;
+    }
+    cout << "output based on "
+         << "sum_ADC = " << sum_ADC << ", zero_sup = "
+         << m_ZeroSuppressionADC << " : ";
+    if (digi_tower)
+    {
+      digi_tower->identify();
+    }
+    else
+    {
+      cout << "None" << endl;
+    }
+  }
+
+  return digi_tower;
+}
+
 void RawTowerDigitizer::CreateNodes(PHCompositeNode *topNode)
 {
   PHNodeIterator iter(topNode);

simulation/g4simulation/g4calo/RawTowerDigitizer.h

@@ -34,15 +34,18 @@ class RawTowerDigitizer : public SubsysReco
   unsigned int get_seed() const { return m_Seed; }
   enum enu_digi_algorithm
   {
-    //! directly pass the energy of sim tower to digitalized tower
+    //! directly pass the energy of sim tower to digitized tower
     kNo_digitization = 0,
     //! wrong spelling, kept for macro compatibility
     kNo_digitalization = 0,
 
-    //! simple digitization with photon statistics, ADC conversion and pedstal
+    //! simple digitization with photon statistics, single amplitude ADC conversion and pedestal
     kSimple_photon_digitization = 1,
     //! wrong spelling, kept for macro compatibility
-    kSimple_photon_digitalization = 1
+    kSimple_photon_digitalization = 1,
+
+    //! digitization with photon statistics on SiPM with an effective pixel N, ADC conversion and pedestal
+    kSiPM_photon_digitization = 2
   };
 
   enu_digi_algorithm
@@ -141,6 +144,12 @@ class RawTowerDigitizer : public SubsysReco
     m_SimTowerNodePrefix = simTowerNodePrefix;
   }
 
+  // ! SiPM effective pixel per tower, only used with kSiPM_photon_digitalization
+  void set_sipm_effective_pixel(const unsigned int &d) { m_SiPMEffectivePixel = d; }
+
+  // ! SiPM effective pixel per tower, only used with kSiPM_photon_digitalization
+  unsigned int get_sipm_effective_pixel() { return m_SiPMEffectivePixel; }
+
  private:
   void CreateNodes(PHCompositeNode *topNode);
 
@@ -149,6 +158,10 @@ class RawTowerDigitizer : public SubsysReco
   //! \return a new RawTower object contain digitalized value of ADC output in RawTower::get_energy()
   RawTower *simple_photon_digitization(RawTower *sim_tower);
 
+  //! digitization with photon statistics on SiPM with an effective pixel N, ADC conversion and pedestal
+  //! this function use the effective pixel to count for the effect that the sipm is not evenly lit
+  RawTower *sipm_photon_digitization(RawTower *sim_tower);
+
   enu_digi_algorithm m_DigiAlgorithm;
 
   RawTowerContainer *m_SimTowers;
@@ -180,6 +193,10 @@ class RawTowerDigitizer : public SubsysReco
   int m_TowerType;
 
   unsigned int m_Seed;
+
+  // ! SiPM effective pixel per tower, only used with kSiPM_photon_digitalization
+  unsigned int m_SiPMEffectivePixel;
+
 #ifndef __CINT__
   gsl_rng *m_RandomGenerator;
 #endif