Changeset 10093 for trunk/Mars/melectronics

Timestamp:

01/07/11 15:48:35 (14 years ago)

Author:

tbretz

Message:

Added afterpulse treatment to MAvalanchePhotoDiode and MSimAPD

File:

• trunk/Mars/melectronics/MAvalanchePhotoDiode.cc (modified) (17 diffs)

trunk/Mars/melectronics/MAvalanchePhotoDiode.cc

@@ -39 +39 @@
 //  Hamamatsu 60x60 cells:  APD(60, 0.2, 3, 8.75)
 //
+//
+// The implementation of afterpulsing is based on
+//    A.Du, F.Retiere
+//    After-pulsing and cross-talk in multi-pixel photon counters
+//    NIM A, Volume 596, Issue 3, p. 396-401
+//
+//
+// Example:
+//
+//    APD apd(ncells, crosstalk, deadtime, recovery);
+//    apd.SetAfterpulseProb(0.14, 0.11);
+//
+//    while (1)
+//    {
+//        // Make the chip "empty" from the influence of external photons
+//        // It also sets fTime to 0.
+//        apd.Init(freq); // freq of external noise (eg. nsb)
+//
+//        // Now call this function for each external photon you have. The
+//        // times are relative to the the time you get by apd.GetTime()
+//        // set automatically after the call to apd.Init().
+//        for (int i=0; i<nphot; i++)
+//            apd.HitRandomCellRelative(dt);
+//
+//        [...]
+//
+//        // Process and produce afterpulses until a time, wrt to GetTime(),
+//        // given by the end of your simulated window. Note that this
+//        // does not produce noise. This must have been properly simulated
+//        // up to this time already.
+//        apd.IncreaseTime(dtend);
+//
+//        // Now you can excess the afterpulses by
+//        TIter Next(&a->GetListOfAfterpulses());
+//        Afterpulse *ap = 0;
+//        while ((ap=static_cast<Afterpulse*>(Next())))
+//        {
+//           if (apd.GetTime()>=dtend)
+//              continue;
+//
+//           cout << "Amplitude:    " << ap->GetAmplitude() << endl;
+//           cout << "Arrival Time: " << ap->GetTime()      << endl;
+//        }
+//     }
+//
+//
 //////////////////////////////////////////////////////////////////////////////
 #include "MAvalanchePhotoDiode.h"
@@ -45 +91 @@
 
 #include "MMath.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
 
 ClassImp(APD);
@@ -77 +126 @@
 APD::APD(Int_t n, Float_t prob, Float_t dt, Float_t rt)
     : fHist("APD", "", n, 0.5, n+0.5, n, 0.5, n+0.5),
-    fCrosstalkProb(prob), fDeadTime(dt), fRecoveryTime(rt), fTime(-1)
+    fCrosstalkProb(prob), fDeadTime(dt), fRecoveryTime(rt),
+    fTime(-1)
 {
     fHist.SetDirectory(0);
-}
-
-// --------------------------------------------------------------------------
-//
-// This is the time the cell needs after a hit until less than threshold
-// (0.001 == one permille) of the signal is lost.
-//
-// If all cells of a G-APD have fired simultaneously and they are fired
-// once more after the relaxation time (with the defaultthreshold of 0.001)
-// the chip will show a signal which is one permille too small.
+
+    fAfterpulses.SetOwner();
+
+    fAfterpulseProb[0] = 0;
+    fAfterpulseProb[1] = 0;
+
+    fAfterpulseTau[0] = 15;
+    fAfterpulseTau[1] = 85;
+}
+
+// --------------------------------------------------------------------------
+//
+// This is the time a chips needs after an external signal to relax to
+// a "virgin" state, i.e. without no influence of the external pulse
+// above the given threshold.
+//
+// It takes into account the dead time of the cell, the relaxation time
+// and the two afterpulse distributions. However, in most cases the
+// afterpulse distribution will dominate (except they are switched off by
+// a zero probability).
+//
+// FIXME: Maybe the calculation of the relaxation time could be optimized?
 //
 Float_t APD::GetRelaxationTime(Float_t threshold) const
 {
-    return fDeadTime-TMath::Log(threshold)*fRecoveryTime;
+    // Calculate time until the probability of one of these
+    // events falls below the threshold probability.
+    const Double_t rt0 = - TMath::Log(threshold)*fRecoveryTime;
+    const Double_t rt1 = fAfterpulseProb[0]>0 ? -TMath::Log(threshold/fAfterpulseProb[0])*fAfterpulseTau[0] : 0;
+    const Double_t rt2 = fAfterpulseProb[0]>0 ? -TMath::Log(threshold/fAfterpulseProb[1])*fAfterpulseTau[1] : 0;
+
+    // Probability not between t and inf, but between t and t+dt
+    // -tau * log ( p / ( 1 - exp(- dt/tau) ) ) = t
+
+    return fDeadTime + TMath::Max(rt0, TMath::Max(rt1, rt2));
 }
 
@@ -108 +179 @@
 // probability.
 //
+// For each photon the possible afterpulses of two distributions are
+// created and added to the list of afterpulses. This is done by calling
+// GenerateAfterpulse for the two afterpulse-distributions.
+//
 // The total height of the signal (in units of photons) is returned.
 // Note, that this can be a fractional number.
@@ -122 +197 @@
     //            y<1 || y>fHist.GetNbinsY())
     //            return 0;
+#ifdef DEBUG
+    cout << "Hit: " << t << endl;
+#endif
 
     // Number of the x/y cell in the one dimensional array
@@ -137 +215 @@
     // Photons within the dead time are just ignored
     if (/*hx.GetBinContent(x,y)>0 &&*/ dt<=0)
+    {
+#ifdef DEBUG
+        cout << "Dead: " << t << " " << cont << " " << dt << endl;
+#endif
         return 0;
-
+    }
     // The signal height (in units of one photon) produced after dead time
     // depends on the recovery of the cell - described by an exponential.
     const Float_t weight = fRecoveryTime<=0 ? 1. : 1-TMath::Exp(-dt/fRecoveryTime);
+
+    // Now we know the charge in the cell and we can generate
+    // the afterpulses with both time-constants
+    GenerateAfterpulse(cell, 0, weight, t);
+    GenerateAfterpulse(cell, 1, weight, t);
 
     // The probability that the cell emits a photon causing crosstalk
@@ -201 +288 @@
 //
 // Check if x and y is a valid cell. If not return 0, otherwise
-// HitCelImp(x, y, t)
+// HitCelImp(x, y, t). HitCellImp generates Crosstalk and Afterpulses.
 //
 // The default time is 0.
@@ -217 +304 @@
 //
 // Determine randomly (uniformly) a cell which was hit. Return
-// HitCellImp for this cell and the given time.
+// HitCellImp for this cell and the given time. HitCellImp
+// generates Crosstalk and Afterpulses
 //
 // The default time is 0.
@@ -238 +326 @@
 // --------------------------------------------------------------------------
 //
-// Sets all cells with a contents whihc is well before the time t such that
+// Sets all cells with a contents which is well before the time t such that
 // the chip is "virgin". Therefore all cells are set to a time which
 // is twice the deadtime before the given time and 1000 times the recovery
 // time.
 //
+// The afterpulse list is deleted.
+//
 // If deadtime and recovery time are 0 then t-1 is set.
 //
@@ -260 +350 @@
     fHist.SetEntries(1);
 
+    fAfterpulses.Delete();
+
     fTime = t;
 }
@@ -267 +359 @@
 // First call FillEmpty for the given time t. Then fill each cell by
 // by calling HitCellImp with time t-gRandom->Exp(n/rate) with n being
-// the total number of cells.
+// the total number of cells. This the time at which the cell was last hit.
 //
 // Sets fTime to t
 //
-// The default time is 0.
+// If the argument t is omitted it defaults to 0.
+//
+// Since after calling this function the chip should reflect the
+// status at the new time fTime=t, all afterpulses are processed
+// until this time. However, the produced random pulses might have produced
+// new new afterpulses.
+//
+// All afterpulses before the new timestamp are deleted.
+//
+// WARNING: Note that this might not correctly reproduce afterpulses
+//          produced by earlier pulese.
 //
 void APD::FillRandom(Float_t rate, Float_t t)
@@ -282 +384 @@
     const Double_t f = (nx*ny)/rate;
 
-    // FIXME: This is not perfect, is it? What about the dead time?
+    // FIXME: Dead time is not taken into account,
+    //        possible earlier afterpulses are not produced.
 
     for (int x=1; x<=nx; x++)
@@ -288 +391 @@
             HitCellImp(x, y, t-MMath::RndmExp(f));
 
+    // Deleting of the afterpulses before fHist.GetMinimum() won't
+    // speed things because we have to loop over them once in any case
+
+    ProcessAfterpulses(fHist.GetMinimum(), t);
+    DeleteAfterpulses(t);
+
     fTime = t;
+}
+
+// --------------------------------------------------------------------------
+//
+// Shift all times including fTime to dt (ie. add -dt to all times)
+// This allows to set a user-defined T0 or shift T0 to fTime=0.
+//
+// However, T0<0 is not allowed (dt cannot be greater than fTime)
+//
+void APD::ShiftTime(Double_t dt)
+{
+    if (dt>fTime)
+    {
+        gLog << warn << "APD::ShiftTime: WARNING - requested time shift results in fTime<0... ignored." << endl;
+        return;
+    }
+
+    // If reset was requested shift all times by end backwards
+    // so that fTime is now 0
+    const Int_t n = (fHist.GetNbinsX()+2)*(fHist.GetNbinsY()+2);
+    for (int i=0; i<n; i++)
+        fHist.GetArray()[i] -= dt;
+
+    fTime -= dt;
 }
 
@@ -296 +429 @@
 // freq. Returns the total signal "recorded".
 //
-// fTime is set to the fTime+dt
+// fTime is set to the fTime+dt.
 //
 // If you want to evolve over a default relaxation time (relax the chip
 // from a signal) use Relax instead.
 //
+// Since after calling this function the chip should reflect the
+// status at the new time fTime=fTime+dt, all afterpulses are processed
+// until this time. However, evolving the chip until this time might
+// have produced new afterpulses.
+//
+// All afterpulses before the new timestamp are deleted.
+//
 Float_t APD::Evolve(Double_t freq, Double_t dt)
 {
-    const Double_t avglen = 1./freq;
-
-    const Double_t end  = fTime+dt;
+    const Double_t end = fTime+dt;
 
     Float_t hits = 0;
-
-    Double_t time = fTime;
-    while (1)
+ 
+    if (freq>0)
     {
-        time += MMath::RndmExp(avglen);
-        if (time>end)
-            break;
-
-        hits += HitRandomCell(time);
+        const Double_t avglen = 1./freq;
+
+        Double_t time = fTime;
+        while (1)
+        {
+            const Double_t deltat = MMath::RndmExp(avglen);
+            if (time+deltat>end)
+                break;
+
+            time += deltat;
+
+            hits += HitRandomCell(time);
+        }
     }
+
+    // Deleting of the afterpulses before fTime won't speed things
+    // because we have to loop over them once in any case
+
+    ProcessAfterpulses(fTime, dt);
+    DeleteAfterpulses(end);
 
     fTime = end;
@@ -329 +480 @@
 // dead.
 // The default time is 0.
+//
+// Note that if you want to get a correct measure of teh number of dead cells
+// at the time t, this function will only produce a valid count if the
+// afterpulses have been processed up to this time.
 //
 Int_t APD::CountDeadCells(Float_t t) const
@@ -348 +503 @@
 // Returs the number of cells which have a time t<=fDeadTime+fRecoveryTime.
 // The default time is 0.
+//
+// Note that if you want to get a correct measure of teh number of dead cells
+// at the time t, this function will only produce a valid count if the
+// afterpulses have been processed up to this time.
 //
 Int_t APD::CountRecoveringCells(Float_t t) const
@@ -364 +523 @@
     return n;
 }
+
+// --------------------------------------------------------------------------
+//
+// Generate an afterpulse originating from the given cell and a pulse with
+// charge. The afterpulse distribution to use is specified by
+// the index. The "current" time to which the afterpulse delay refers must
+// be given by t.
+//
+// A generated Afterpulse is added to the list of afterpulses
+//
+void APD::GenerateAfterpulse(UInt_t cell, Int_t idx, Double_t charge, Double_t t)
+{
+    // The cell had a single avalanche with signal height weight.
+    // This cell now can produce an afterpulse photon/avalanche.
+    const Double_t p = gRandom->Uniform();
+
+    // It's probability scales with the charge of the pulse
+    if (p>charge*fAfterpulseProb[idx])
+        return;
+
+    // Afterpulses come with a well defined time-constant
+    // after the normal pulse
+    const Double_t dt = MMath::RndmExp(fAfterpulseTau[idx]);
+
+    fAfterpulses.Add(new Afterpulse(cell, t+dt));
+
+#ifdef DEBUG
+    cout << "Add : " << t << " + " << dt << " = " << t+dt << endl;
+#endif
+}
+
+// --------------------------------------------------------------------------
+//
+// Process afterpulses between time and time+dt. All afterpulses in the list
+// before t=time are ignored. All afterpulses between t=time and
+// t=time+dt are processed through HitCellImp. Afterpulses after and
+// equal t=time+dt are skipped.
+//
+// Since the afterpulse list is a sorted list newly generated afterpulses
+// are always inserted into the list behind the current entry. Consequently,
+// afterpulses generated by afterpulses will also be processed correctly.
+//
+// Afterpulses with zero amplitude are deleted from the list. All other after
+// pulses remain in the list for later evaluation.
+//
+void APD::ProcessAfterpulses(Float_t time, Float_t dt)
+{
+#ifdef DEBUG
+    cout << "Process afterpulses from " << time << " to " << time+dt << endl;
+#endif
+
+    const Float_t end = time+dt;
+
+    TObjLink *lnk = fAfterpulses.FirstLink();
+    while (lnk)
+    {
+        Afterpulse &ap = *static_cast<Afterpulse*>(lnk->GetObject());
+
+        // Skip afterpulses which have been processed already
+        // or which we do not have to process anymore
+        if (ap.GetTime()<time || ap.GetAmplitude()>0)
+        {
+            lnk = lnk->Next();
+            continue;
+        }
+
+        // No afterpulses left in correct time window
+        if (ap.GetTime()>=end)
+           break;
+
+        // Process afterpulse through HitCellImp
+        const Float_t ampl = ap.Process(*this);
+
+        // Step to the next entry already, the current one might get deleted
+        lnk = lnk->Next();
+
+        if (ampl!=0)
+            continue;
+
+#ifdef DEBUG
+        cout << "Del : " << ap.GetTime() << " (" << ampl << ")" << endl;
+#endif
+
+        // The afterpulse "took place" within the dead time of the
+        // pixel ==> No afterpulse, no crosstalk.
+        delete fAfterpulses.Remove(&ap);
+    }
+}
+
+// --------------------------------------------------------------------------
+//
+// Delete all afterpulses before and equal to t=time from the list
+//
+void APD::DeleteAfterpulses(Float_t time)
+{
+    TIter Next(&fAfterpulses);
+
+    Afterpulse *ap = 0;
+
+    while ((ap = static_cast<Afterpulse*>(Next())))
+    {
+        if (ap->GetTime()>=time)
+            break;
+
+        delete fAfterpulses.Remove(ap);
+    }
+}