Changeset 10093

Timestamp:

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

Author:

tbretz

Message:

Added afterpulse treatment to MAvalanchePhotoDiode and MSimAPD

Location:

trunk/Mars

Files:

• Changelog (modified) (1 diff)
• melectronics/MAvalanchePhotoDiode.cc (modified) (17 diffs)
• msimcamera/MSimAPD.cc (modified) (9 diffs)

trunk/Mars/Changelog

@@ -18 +18 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+ 2011/01/07 Thomas Bretz
+
+   * msimcamera/MSimAPD.[h,cc]:
+     - added setting of the afterpulse probabilities
+     - handle afterpulses in processing
+
+   * melectronics/MAvalanchePhotoDiode.[h,cc]:
+     - restructured header file
+     - added treatment of afterpulses in HitCellImp
+     - added time constants and afterpulse probability 
+       data memebers
+     - added TSortedList to hold the afterpulses
+     - added afterpulses to calculation of relaxation time
+     - added call to Relax to Init
+     - added treatment of afterpulses to FillEmpty, FillRandom and
+       Evolve
+     - added member functions to process the afterpulses
+
+
 
  2011/01/06 Thomas Bretz

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);
+    }
+}

trunk/Mars/msimcamera/MSimAPD.cc

@@ -37 +37 @@
 //
 // Remark:
-//   - The photons MUST be sorted increasing in time.
 //   - The photon rate used to initialize the APD must match the one used
 //     to "fill" the random photons. (FIXME: This should be stored somewhere)
@@ -152 +151 @@
     //   * Make the arguments a data member of MSimAPD
 
-    Int_t   ncells    = 0;
-    Float_t crosstalk = 0;
-    Float_t deadtime  = 0;
-    Float_t recovery  = 0;
+    Int_t   ncells     = 0;
+    Float_t crosstalk  = 0;
+    Float_t deadtime   = 0;
+    Float_t recovery   = 0;
+    Float_t afterprob1 = 0;
+    Float_t afterprob2 = 0;
 
     switch (fType)
     {
     case 1:
-        ncells    = 30;
-        crosstalk = 0.2;
-        deadtime  = 3;
-        recovery  = 8.75*4;
+        ncells     = 30;
+        crosstalk  = 0.2;
+        deadtime   = 3;
+        recovery   = 8.75*4;
+        afterprob1 = 0;
+        afterprob2 = 0;
         break;
 
     case 2:
-        ncells    = 60;
-        crosstalk = 0.2;
-        deadtime  = 3;
-        recovery  = 8.75;
+        ncells     = 60;
+        crosstalk  = 0.2;
+        deadtime   = 3;
+        recovery   = 8.75;
+        afterprob1 = 0;
+        afterprob2 = 0;
         break;
 
     case 3:
-        ncells    = 60;
-        crosstalk = 0.15;
-        deadtime  = 3;
-        recovery  = 8.75;
+        ncells     = 60;
+        crosstalk  = 0.15;
+        deadtime   = 3;
+        recovery   = 8.75;
+        afterprob1 = 0;
+        afterprob2 = 0;
+        break;
+
+    case 4:
+        ncells     = 60;
+        crosstalk  = 0.15;
+        deadtime   = 3;
+        recovery   = 8.75;
+        afterprob1 = 0.14;
+        afterprob2 = 0.11;
         break;
 
@@ -186 +202 @@
 
     for (UInt_t i=0; i<fGeom->GetNumPixels(); i++)
-        fAPDs.Add(new APD(ncells, crosstalk, deadtime, recovery));
+    {
+        APD *apd = new APD(ncells, crosstalk, deadtime, recovery);
+        apd->SetAfterpulseProb(afterprob1, afterprob2);
+
+        fAPDs.Add(apd);
+    }
 
     return kTRUE;
@@ -198 +219 @@
 Int_t MSimAPD::Process()
 {
-    //const Double_t rate = 40e9;
-    // FIXME: Where do we get this number from??
-    // const Double_t rate = 0.04;
-
     // Make all APDs look neutral for the first hit by a photon according to the
     // average hit rate
@@ -213 +230 @@
         return kERROR;
     }
-/*
-    for (UInt_t idx=0; idx<npix; idx++)
-    {
-        const Double_t freq = fRates ? (*fRates)[idx].GetPedestal() : fFreq;
-        static_cast<APD*>(fAPDs.UncheckedAt(idx))->FillRandom(freq, fStat->GetTimeFirst());
-    }
-*/
+
+    // To ensure that the photons are always sorted when calling
+    // HitRandomCellRelative the array is sorted here. If it is sorted
+    // already nothing will be done since the status is stored.
+    // FIXME: Check that this is true and check that it is really necessary
+    fEvt->Sort();
 
     // This tries to initialize dead and relaxing cells properly. If
@@ -228 +244 @@
     {
         const Double_t freq = fRates ? (*fRates)[idx].GetPedestal() : fFreq;
+
+        // Init creates an empty G-APD, i.e. without external pulses
+        // but the correct history for afterpulses and relaxation.
+        // If it was already initialized once it evolves the G-APD
+        // for a time until the effect of relaxation and afterpulses
+        // is below 0.1%. The also creates the possible afterpulses
+        // of the future and deletes later afterpulses from the list.
+        // After the the time stamp fTime is set to 0.
         static_cast<APD*>(fAPDs.UncheckedAt(idx))->Init(freq);
     }
@@ -240 +264 @@
         MPhotonData &ph = (*fEvt)[i];
 
-        // Get arrival time of photon and idx
+        // Get arrival time of photon wrt to left edge of window and its index
         const Double_t t = ph.GetTime()-fStat->GetTimeFirst();
         const Int_t  idx = ph.GetTag();
@@ -255 +279 @@
             return kERROR;
         }
-        // Simulate hitting the APD (the signal height in effective
-        // "number of photons" is returned)
+
+        // Simulate hitting the APD at a time t after T0 (APD::fTime).
+        // Crosstalk is taken into account and the resulting signal height
+        // in effective "number of photons" is returned. Afterpulses until
+        // this time "hit" the G-APD and newly created afterpulses
+        // are stored in the list of afterpulses
         const Double_t hits = static_cast<APD*>(fAPDs.UncheckedAt(idx))->HitRandomCellRelative(t);
-
-        // FIXME: Make a proper simulation of the excess noise!!!
-        //const Double_t signal = gRandom->Gaus(hits, 0.2*TMath::Sqrt(hits));
 
         // Set the weight to the input
@@ -269 +294 @@
     // simulated time.
     for (UInt_t idx=0; idx<npix; idx++)
-        static_cast<APD*>(fAPDs.UncheckedAt(idx))->IncreaseTime(fStat->GetTimeLast());
+    {
+        APD *a = static_cast<APD*>(fAPDs.UncheckedAt(idx));
+
+        const Double_t end = fStat->GetTimeLast()-fStat->GetTimeFirst();
+
+        // This moves T0 (APD::fTime) at the right edge of our time-
+        // window. For the next event this means that afterpulses of past
+        // noise and afterpulses will be available already.
+        // FIXME: Note, that this might mean that a cosmics-pulse
+        //        might increase the noise above the normal level.
+        a->IncreaseTime(end);
+
+        // Get the afterpulses and add them to the signal
+        TIter Next(&a->GetListOfAfterpulses());
+        Afterpulse *ap = 0;
+        while ((ap=static_cast<Afterpulse*>(Next())))
+        {
+            // Skip afterpulses later than that which have been
+            // already produced
+            if (ap->GetTime()>=end)
+                continue;
+
+            // Add a new photon
+            // FIXME: SLOW!
+            MPhotonData &ph = fEvt->Add();
+
+            // Set source to Artificial (noise), the amplitude produced by the
+            // afterpulse (includes crosstalk), its arrival time
+            // and its amplitude, as well as the index of the channel it
+            // corresponds to.
+            ph.SetPrimary(MMcEvtBasic::kArtificial);
+            ph.SetWeight(ap->GetAmplitude());
+            ph.SetTime(ap->GetTime()+fStat->GetTimeFirst());
+            ph.SetTag(idx);
+        }
+
+        // It seems to make sense to delete the previous afterpulses now
+        // but this is not necessary. We have to loop over them in any
+        // case. So we omit the additional loop for deletion but instead
+        // do the deletion in the next loop at the end of Init()
+        // If needed this could be used to keep the total memory
+        // consumption slighly lower.
+    }
+
+    // Now the newly added afterpulses have to be sorted into the array correctly
+    fEvt->Sort();
 
     return kTRUE;