scintillation time function changed

include/utils.h

@@ -162,12 +162,12 @@ double AttenuationFactor(double d, int planeID);
 
 namespace scintillation
 {
- const double tscin = 3.08;
- const double tscex = 0.588;
- const double vlfb = 5.85;
- const double rise_time = 0.7; // ns
- const double decay_time = 3.0; // ns
-} 
+const double tscin = 3.08;
+const double tscex = 0.588;
+const double vlfb = 5.85;
+const double rise_time = 0.7; // ns
+const double decay_time = 3.0; // ns
+}  // namespace scintillation
 
 namespace photo_sensor
 {

src/SANDDigitization.cpp

@@ -13,30 +13,49 @@ TRandom3 rand(0);
 
 namespace ecal
 {
+
+// simulate time of emission of the scintillation photons
+double photons_emission_time()
+{
+ // https://github.com/Geant4/geant4/blob/e58e650b32b961c8093f3dd6a2c3bc917b2552be/source/processes/electromagnetic/xrays/src/G4Scintillation.cc#L638
+ double t;
+ do {
+ // The exponential distribution as an envelope function: very efficient
+ t = -1.0 * sand_reco::ecal::scintillation::decay_time * TMath::Log(1.0 - rand.Uniform() );
+ } while (rand.Uniform() >
+ (1.0 - TMath::Exp(-t / sand_reco::ecal::scintillation::rise_time)));
+ return t;
+}
+
+double photons_emission_time_as_kloe()
+{
+ return sand_reco::ecal::scintillation::tscin *
+ TMath::Power(1. / rand.Uniform() - 1.,
+ sand_reco::ecal::scintillation::tscex);
+}
+
 // simulate pe arrival time to pmt
 double photo_electron_time_to_pmt_arrival_time(double t0, double d)
 {
  /*
- - For each photoelectron: Time for TDC simulation obtained from
-
-C PHOTOELECTRON TIME : Particle TIME in the cell
+ For each photoelectron: Time for TDC simulation obtained from
+ C PHOTOELECTRON TIME : Particle TIME in the cell
 C + SCINTILLATION DECAY TIME +
 C + signal propagation to the cell
 C + 1ns uncertainty
-
- TPHE = Part_time+TSDEC+DPM1*VLFB+Gauss(1ns)
-
- VLFB = 5.85 ns/m
+ TPHE = Part_time+TSDEC+DPM1*VLFB+Gauss(1ns)
+VLFB = 5.85 ns/m
 !!!! Input-TDC Scintillation time -
- TSDEC = TSCIN*(1./RNDMPH(1)-1)**TSCEX (ns)
+ TSDEC = TSCIN*(1./RNDMPH(1)-1)**TSCEX (ns)
+TSCIN 3.08 ns
+TSCEX 0.588
+*/
 
-  TSCIN 3.08 ns
- TSCEX 0.588
- */
+ // double tdec = sand_reco::ecal::scintillation::tscin *
+ //  TMath::Power(1. / rand.Uniform() - 1.,
+ // sand_reco::ecal::scintillation::tscex);
 
- double tdec = sand_reco::ecal::scintillation::tscin *
- TMath::Power(1. / rand.Uniform() - 1.,
- sand_reco::ecal::scintillation::tscex);
+ double tdec = photons_emission_time();
 
  double time = t0 + tdec +
  sand_reco::ecal::scintillation::vlfb * d * conversion::mm_to_m +